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I Introduction 


The  implementation  of  modern  control  techniques  re- 
quires a detailed  knowledge  of  the  controlled  system.  To 
accomplish  this  objective,  a mathematical  model  of  the  system 
is  developed,  but  this  model  generally  contains  parameters 
which  are  known  only  with  some  limited  degree  of  certainty. 

It  may  be  necessary  for  adequate  control  that  these  unknown 
parameters  be  estimated  as  accurately  as  possible. 

One  method  of  obtaining  the  estimates  is  for  the  engi- 
neer, based  on  his  experience  and  familiarity  with  the  sys- 
tem, to  assign  a value  to  the  parameters.  This  is  usually 
not  adequate  for  most  system  functions.  A better  method  is 
to  perform  experiments  on  the  system  and,  by  observing  the 
inputs,  outputs,  and  other  available  data,  to  estimate  the 
parameters.  It  is  well  known  that  the  choice  of  the  input 
signals  has  a direct  bearing  on  the  accuracy  of  the  param- 
eter estimates  for  most  systems.  The  objective  of  this  re- 
search is  to  improve  the  estimates  of  the  parameters  by 
modifying  the  input  controls. 

A typical  controlled  system  is  shown  in  Figure  1. 
Uncertainties  and  disturbances  are  also  generally  inputs  to 
the  system,  but  are  not  addressed  directly  in  this  research. 
The  system  output  measurements  are  corrupted  by  noise  and 
the  input  control  signals  usually  belong  to  either  an  open- 
loop,  feedback,  or  closed-loop  policy.  The  usual  definition 
of  these  different  controls  as  given  by  Bar-Shalom  and  Tse 
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Noise 


Input 

Controls 


■*> 


Output 

Measurements 


Figure  1 . Controlled  System  Diagram 

(Ref  5)  will  be  used  in  this  research.  An  open-loop  con- 
trol policy  is  determined  from  a priori  information  only. 
Even  if  output  measurements  are  taken,  this  information  is 
not  used  to  calculate  the  controls.  The  complete  control 
sequence  is  determined  before  the  test  has  begun.  A feed- 
back control  policy  is  calculated  from  a priori  information 
plus  any  measured  data  up  to  and  including  the  time  the 
control  is  required.  It  ignores  the  fact  that  future  mea- 
surements will  be  made,  while  a closed-loop  control  policy 
uses  the  same  information  as  the  feedback  plus  it  exploits 
the  fact  that  future  measurements  will  be  made  in  its  cal- 
culation of  appropriate  controls. 

It  is  shown  in  the  next  section  that  most  research  has 
addressed  calculating  an  optimal  open-loop  control  policy 
that  improves  the  identification  of  unknown  parameters  in 
the  system.  The  advantage  of  the  open-loop  controls  is 
that  they  are  usually  easier  to  calculate  than  feedback  and 
closed-loop  controls.  An  expected  disadvantage  is  that 
noise  and  other  uncertainties  modify  the  system's  response 
to  the  controls,  so  that  the  control  sequence  is  no  longer 


optimal. 

A method  of  monitoring  the  output  and  adjusting  the  in- 
put control  accordingly  is  provided  by  the  use  of  a feedback 
term.  The  advantages  of  feedback  to  control  the  output  with 
parameter  changes  is  well  known,  but  little  work  has  been 
done  in  using  feedback  to  enhance  the  information  in  the  out- 
put about  parameter  changes.  This  dissertation  develops  a 
method  of  computing  a feedback  control  policy  that  improves 
the  identification  of  unknown  parameters  in  the  system  beyond 
that  achievable  with  purely  open-loop  inputs.  The  control 
will  be  of  the  discrete-time  form 

u(k)  = Fyi(k)  ♦ v(k)  (1) 

where  u(k)  is  an  m-dimensional  input  control  at  the  k-th 
time  instant,  F is  a constant  mxr  feedback  matrix,  ^(k)  is 
the  r-dimensional  output  vector  at  the  k-th  sample  time  in- 
stant, and  v(k)  is  an  m-dimensional  open-loop  control  vector 
at  the  k-th  instant.  The  terms  Fy  and  v(k)  will  be  calcu- 
lated a priori;  however,  u(k)  is  a feedback  control  policy 
because  of  the  output  vector,  £.(k),  appearing  in  Eq  (1). 

Figure  2.  shows  the  modified  controlled  system.  The 
addition  of  this  feedback  terra  increases  the  degree  of 
freedom  we  have  to  optimize  the  information  available  in 
the  output  measurements  for  parameter  identification.  It 
is  a means  of  adjusting  the  system  eigenvalues.  This  may 
be  necessary  in  order  to  stabilize  an  originally  unstable 
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Figure  2.  Controlled  System  with  Feedback 

system  or  to  make  an  originally  unidentifiable  system  iden- 
tifiable. It  is  shown  in  this  dissertation  that  the  addi- 
tion of  feedback  of  the  output  measurements  greatly  improves 
the  parameter  estimates  when  in  the  presence  of  a high  mea- 
surement noise  environment. 

To  limit  the  scope  of  the  dissertation,  all  parameters 
are  assumed  constant,  and  a nominal  value  for  each  parameter 
is  known  a priori.  Only  linear  time-invariant  mathematical 
models  of  systems  are  addressed.  These  models  are  the  para- 
metric state  models,  and  the  initial  conditions  on  the  state 
vector  is  assumed  to  be  either  zero  or  calculated  by  the  pro- 
cedure developed  in  the  next  chapter.  This  implies  that 
there  is  an  optimal  initial  condition  for  the  state  vector 
to  improve  the  parameter  estimates.  Dynamic  driving  noise 


is  not  considered  in  this  research.  It  will  be  shown  in  a 
later  section  of  this  chapter  that  the  addition  of  process 
noise  increases  the  complexity  of  the  problem. 

Background  Literature 

A considerable  amount  of  research  has  been  performed  in 
the  area  of  optimal  inputs  for  parameter  identification.  The 
word  "optimal"  implies  that  some  criterion  must  be  used  to 
select  the  desired  control.  The  criterion  used  in  most  of 
the  research  is  either  maximization  of  a scalar  function  of 
the  Fisher  information  matrix  (Ref  27)  or  minimization  of  a 
scalar  function  of  the  parameter  error  covariance  matrix. 

The  use  of  the  parameter  error  covariance  matrix  seems 
natural  since  the  primary  purpose  is  to  obtain  the  most  ac- 
curate parameter  estimates.  Using  controls  to  minimize 
functions  of  this  error  matrix  should  give  us  these  accu- 
rate estimates.  Unfortunately,  in  general,  it  is  mathemat- 
ically complicated  to  use  this  criterion;  therefore,  the 
use  of  the  Fisher  information  matrix  became  attractive. 

The  information  matrix  is  a measure  of  the  amount  of 
information  about  the  parameters  tnat  is  available  in  the 
output.  Maximizing  this  matrix  in  some  manner  is  then  a 
means  of  increasing  the  information  available  about  the  un- 
known parameters.  However,  "maximizing  a matrix"  usually 
is  accomplished  by  maximizing  a scalar  function  of  the  ma- 
trix. Some  of  the  function  used  have  been  the  trace,  the 
weighted  trace,  and  the  determinant  of  the  matrix.  Some 
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advantages  and  disadvantages  of  using  these  various  func- 
tions will  be  presented  when  they  are  first  discussed  in  the 
following  paragraphs. 
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Most  of  the  previous  research  has  been  in  the  area  of 
developing  open-loop  inputs  for  both  linear  and  nonlinear 
systems  with  measurement  noise.  Process  noise  is  usually 
not  considered.  The  main  reasons  are  that  it  results  in  a 
more  difficult  math  model  of  the  system  and  the  computation 
required  is  more  involved,  '//hen  process  noise  is  used,  the 
system  considered  is  usually  linear  and  described  by  a low 
dimension  model.  Levin  (Ref  16)  considers  single-input, 
single-output  linear  systems  without  process  noise.  His 
criterion  is  based  on  the  parameter  error  covariance  matrix. 
The  parameters  that  are  estimated  represent  the  impulse  re- 
sponse of  the  system.  His  system  model  is,  therefore,  lin- 
ear in  the  parameters  so  he  uses  the  least  square  estimation 
technique.  He  shows  that  the  parameter  error  covariance  ma- 
trix is  made  up  of  input  control  terms  only  and  that  the  in- 
verse of  this  matrix  has  identical  terms  on  the  main  diago- 
nal. Levin  uses  a principal,  that  he  proves  in  the  article, 
that  by  making  all  the  off  diagonal  terms  of  the  inverse  ma- 
trix zero,  he  minimizes  the  diagonal  terms  of  the  error  ma- 
trix. For  a stationary  random  process  input  with  zero  mean, 
he  shows  that  this  implies  that  the  input  must  be  white. 

Mehra  (Ref  20  and  21)  considers  linear  and  nonlinear 
systems  with  and  without  process  noise.  His  criterion  is 
the  maximizing  of  the  weighted  trace  of  the  Fisher  informa- 


tion  matrix.  The  advantage  of  using  the  trace  is  that  it  is 
easier  to  compute  than  the  determinant,  but  the  main  disad- 
vantage is  that  this  criterion  could  lead  to  cases  in  which 
the  system  is  unidentifiable.  This  along  with  a discussion 
of  the  possible  criteria  is  given  in  the  section  on  multiple 
unknown  parameters  in  Chapter  IV.  To  overcome  this  problem, 
the  information  matrix  is  multiplied  by  a weighting  matrix 
in  which  the  weighting  factors  are  properly  selected  to  keep 
the  system  identifiable.  Mehra  (Ref  21)  and  Gupta  and  Hall 
(Ref  11)  have  shown  how  to  select  these  weighting  factors  in 
such  a way  that  the  resulting  optimal  control  minimizes  ei- 
ther the  determinant  or  trace  of  the  inverse  of  the  informa- 
tion matrix. 

Aoki  and  Staley  (Ref  1)  basically  consider  the  same 
problem.  Mehra  considers  energy  constraints  on  the  inputs 
and  similarly,  Aoki  and  Staley  also  consider  energy  con- 
straints on  the  output. 

Nahi  and  Wallis  (Ref  23),  Wallis  (Ref  34),  and  Nap jus 
(Ref  24)  consider  nonlinear  systems  without  process  noise. 
Their  optimality  criterion  is  the  maximizing  of  the  weighted 
trace  of  the  information  matrix.  To  calculate  the  informa- 
tion matrix  they  use  augmented  sensitivity  equations,  which 
are  equations  relating  the  sensitivity  of  the  output  to  pa- 
rameter changes.  Amplitude  constraints  on  the  inputs  are 
also  incorporated. 

Reid  (Ref  28)  and  Goodwin  and  Payne  (Ref  9)  use  the 
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minimization  of  the  weighted  inverse  of  the  information  ma- 
trix as  their  criterion.  Using  this  criterion  is  more  dif- 
ficult than  maximizing  the  trace  of  the  information  matrix, 
but  it  does  not  lead  to  the  problem  of  unidentifiable  param- 
eters. 

All  of  the  work  referenced  above  produce  open-loop  con- 
trols that  are  based  on  a priori  information  only.  It  is 
found  that  in  cases  where  the  output  noise  is  large  or  the 
initial  estimates  of  the  unknown  parameters  are  very  erro- 
neous,  the  identification  process  produces  very  inaccurate 
estimates.  Gustavsson,  Ljung,  and  Soderstrora  state  in  their 
survey  paper  (Ref  12)  that  feedback  terms  may  have  to  be  in- 
cluded in  the  input  signal,  feedback  plus  open-loop  controls, 
in  order  to  obtain  acceptable  accuracy.  For  large  distur- 
bances and  high  noise  levels,  feedback  terms  may  be  needed. 
Inclusion  of  feedback  will  also  increase  the  degree  of  free- 
dom when  choosing  the  input  signals.  They  further  state 
that  it  is  difficult  to  determine  a priori  whether  a feed- 
back term  is  necessary  to  give  better  estimates;  therefore, 
the  optimal  inputs  should  include  feedback  terms  with  the 
open-loop  terms  instead  of  open-loop  controls  only.  Al- 
though the  addition  of  feedback  potentially  enhances  the 
parameter  identification  significantly,  the  complexity  of 
solving  this  generalized  problem  is  substantially  greater 
than  that  associated  with  open-loop  controls  only. 

In  an  attempt  to  develop  a closed-loop  control,  Keviczky 
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and  Banyasz  (Ref  15)  and  Arimoto  and  Kimura  (Ref  3)  approach 
the  identification  problem  by  designing  inputs  that  maximize 
the  incremental  increase  in  information  during  the  next  mea- 
surement period.  The  inputs  are  not  true  closed-loop  con- 
trols since  they  do  not  take  future  learning  into  account. 
The  one  step  ahead  approach  is  used  to  simplify  the  calcu- 
lations so  that  an  on-line  algorithm  could  be  used.  Both 
articles  consider  an  amplitude  constraint  on  the  input  sig- 
nal. 

Lopez-Toledo  (Ref  17)  attempts  to  find  the  optimal 
closed-loop  input  with  an  energy  constraint  that  will  max- 
imize the  trace  of  the  expectation  of  the  inverse  of  the 
conditional  error  covariance  matrix.  He  then  showed  that 
for  his  model,  linear  in  the  parameters,  the  "curse  of  di- 
mensionality" makes  the  problem  intractable.  He  then  de- 
velops suboptimal  inputs  by  using  an  affine  law  (input  is 
a linear  function  of  the  output  plus  a constant  term)  and 
then  the  open-loop  feedback  optimal  (OLFO)  approach.  The 
OLFO  method  is  to  compute  an  optimal  open-loop  control  at 
each  time  step,  but  apply  only  the  first  control  and  then 
recompute  at  the  next  time  step.  Lope z-Tole do's  affine  law 
is  the  same  as  the  control  policy  used  in  this  research. 

The  main  differences  are  that  this  research  does  not  re- 
strict the  model  to  be  linear  in  the  parameters,  but  does 
restrict  the  freedom  of  the  closed-loop  system  eigenvalues. 
These  differences  were  chosen  because  most  real  systems  can- 


not  be  put  into  the  linear-in-parameters  form  and  adding 
feedback  moves  the  eigenvalue!  which  must  be  constrained  to 
maintain  stability. 


Dpadhyaya  and  Sorenson  (Ref  32  and  33)  develop  a feed- 
back control  by  casting  the  problem  as  an  OLFO  control  prob- 
lem. The  model  used  is  linear  in  the  parameters,  and  when 
the  a priori  distribution  of  the  parameters  is  known,  the 
optimal  signal  is  obtained  by  maximizing  the  Bhattacharyya 
coefficients,  3^,  where  i ? j.  The  coefficients  are  de- 
fined as 

3i.  = -lnjr(p(zk|ei)p(zk|e.))^dzk  (2) 

id 

where  p(Z  | 9^)  is  the  probability  distribution  of  the  mea- 
surement  sequence,  Z , given  the  parameter  vector,  0^.  . 0 

is  assumed  to  take  on  only  a discrete  number  of  values  so 
Bfj  is  finite  dimensional.  Maximizing  the  Bhattacharyya 
measure  maximizes  the  distance  between  the  pairs  of  dis- 
tributions and  makes  the  identification  process  more  accu- 
rate. The  resultant  control  is  a feedback  control. 

Goodwin  and  Payne  (Ref  10)  develop  optimal  open-loop 
and  feedback  controls  by  minimizing  the  negative  of  the  log 
of  the  determinant  of  the  information  matrix.  They  show 
that  if  there  are  no  unknown  parameters  that  are  in  both 
the  plant  matrix  and  measurement  noise  matrix  and  the  total 
input  power,  open-loop  and  feedback  controls,  is  constrained, 
then  the  optimal  control  consists  of  only  an  open-loop  con- 
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trol.  In  this  research,  however,  the  energy  constraint  is 
only  on  the  open-loop  controls  and  not  on  the  combined  con- 
trols. 


The  next  section  discusses  the  reasoning  behind  selec- 
ting a control  made  up  of  the  sum  of  a linear  feedback  terra 
and  an  open-loop  term.  The  section  also  investigates  the 
different  models  available  and  discusses  why  the  linear, 
time-invariant  model  with  measurement  noise  only  was  selec- 
ted for  this  research. 

Problem  Formulation 

Consider  a linear  time-invariant  discrete  system 


x(k+1)  = A(a)x(k)  + B(a)u(k) 

(3) 

£(k)  = Cx(k) 

U) 

where  x(k)  is  an  n-dimensional  state  vector  at  time  instant 
k,  u(k)  is  an  m-dimensional  control  vector,  and  ^(k)  is  an 
r-dimensional  output  vector.  A(a)  is  an  nxn  time-invariant 
plant  matrix,  B(a)  is  an  nxm  time-invariant  control  matrix, 
and  C is  an  rxn  time-invariant  output  matrix.  The  symbol  a 
represents  the  p-dimensional  unknown  parameter  vector.  As 
mentioned  earlier,  it  is  assumed  that  the  parameters  are 
constant  and  that  a nominal  value  of  the  parameter  vector  a 
is  known  and  is  represented  by  aQ.  The  plant  and  control 
matrices  are  a function  of  a while  the  output  matrix  is  not. 
This  is  a reasonable  constraint  since  the  uncertainties  in 
sensors  ore  usually  less  than  in  a dynamic  system  or  the 


actuators.  The  matrices  are  not  normally  linear  with  re- 
spect to  the  unknown  parameters  when  they  represent  a phys- 
ical system.  This  also  occurs  when  a continuous  system  is 
changed  to  a discrete  system  and  unknown  parameters  exist  in 
the  plant  matrix  of  the  continuous  system.  In  the  discrete 
system  the  unknown  parameter  will  usually  appear  nonlinearly 
in  the  plant  and  control  matrices.  The  system  in  this  re- 
search is  considered  to  be  both  completely  observable  and 
completely  controllable.  The  initial  condition  of  the  state 
vector,  x(0),  will  be  known.  Its  value  will  be  either  0 or 
will  be  calculated  to  optimize  the  estimation  technique. 

The  control  is  chosen  to  have  the  form 


u(k)  = Fy£(k)  + v(k)  (1) 

that  is,  the  control  is  a combination  of  a linear  function 
of  the  current  output  and  an  open-loop  control.  If  F = 0, 
then  the  optimal  input  control  sequence  u(k)  could  be  cal- 
culated from  some  of  the  referenced  methods  for  open-loop 
controls.  If  v(k)  = 0,  then  only  the  feedback  term  is  left. 
Using  the  feedback  control  u(k)  = F £(k)  to  identify  the 
system  could  give  erroneous  results.  Box  and  MacGregor 
(Ref  6)  show  that  identification  techniques  that  are  based 
on  this  input  control,  u(k),  and  the  output,  ^(k) , may  iden- 
tify the  pseudoinverse  of  the  feedback  matrix,  that  is,  F *. 

~y 

Wellstead  (Ref  35)  shows  that  the  difficulties  associated 
with  identification  of  systems  with  feedback  controls  are 
most  readily  solved  by  the  addition  of  an  external  signal. 


Gustavsson,  Ljung,  and  Soderstrom  (Ref  12)  state  that  if  the 
number  of  independent  "extra  inputs",  v(k),  is  equal  to  the 
number  of  inputs  to  the  system,  then  the  system  is  always 
identifiable  no  matter  what  the  feedback  is.  It  is  for 
these  reasons  that  the  input  control  used  in  this  research 
is  of  the  form  of  Eq  (1). 

The  feedback  matrix  will  move  the  closed-loop  system 

poles,  the  eigenvalues  of  A(a)  + B(a)F,C  , to  enhance  the 

y- 

identification.  Movement  of  the  poles  also  indicates  the 
possibility  of  causing  the  feedback  system  to  be  unstable. 
Constraints  must  be  placed  on  these  poles  in  order  to  pre- 
vent  this  from  happening.  The  open-loop  control,  v(k),  then 
not  only  optimizes  the  identification,  but  also  eliminates 
the  danger  of  unidentifiable  parameters  as  mentioned  in  this 
section. 

The  criterion  that  is  used  is  the  maximization  of  a 
function  of  the  Fisher  information  matrix.  This  criterion 
is  chosen  because  as  mentioned  in  the  previous  section,  it 
is  mathematically  easier  to  compute  and  maximizes  the  infor- 
mation about  the  parameters  in  the  output.  The  information 
matrix  is 

" 3 In  p(YJa)T  3 In  p( Y„ j S)~ 

M = E -A  ~ (5) 

YN  L 3 | 3 I 

where  is  a set  of  all  measurements  up  to  time  instant  N, 
p(b|c)  is  the  probability  distribution  of  random  variable  b 
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given  that  c has  occured,  and  E | b represents  the  expec- 

Y L J 

tation  of  b taken  over  the  sample  space  of  observations  Y^. 
Appendix  B outlines  the  convention  that  is  used  when  the 
partial  of  a scalar  is  taken  with  respect  to  a vector.  The 
conditional  likelihood  function  In  pCY^a)  must  be  evaluated, 
and  it  is  this  function  that  has  a bearing  on  the  system 
model  selected  in  this  research. 

The  function  of  the  information  matrix  that  is  maxi- 
mized in  this  research  is  the  trace  or  'weighted  trace.  A 
weighted  trace  means  that  the  information  matrix  is  premul- 
tiplied by  a matrix  whose  elements  are  weighting  factors  and 
then  the  trace  operation  is  performed.  These  weighting  fac- 
tors can  be  determined  a priori  or  during  the  iterative  pro- 
cess of  computing  the  optimal  controls.  Experience  may  in- 
dicate that  information  in  certain  outputs  is  more  important 
than  in  others,  so  these  terms  'would  receive  a higher  weight- 
ing factor  than  the  others.  The  engineer  could  assign  these 
values  a priori;  Chapter  IV  gives  a method  for  computing  the 
factors  during  the  iteration  steps,  but  as  will  be  shown, 
this  requires  much  more  computation  time. 

Equations  (3)  and  (4)  are  deterministic  equations  that 
are  assumed  to  model  the  system.  With  no  process  noise  and 
x(0)  known,  the  output  vector  ^(k)  can  be  measured  exactly 
for  any  input  control  sequence  u(k-1),  u(k-2),***  , u(0). 

The  only  unknown  in  the  equations  is  a.  If  the  control  se- 
quence is  such  that  the  system  is  identifiable  (that  is,  a 
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can  be  uniquely  determined  from  a given  set  of  experimental 
observations)  then  a can  be  found  by  a nonlinear  optimiza- 
tion method  such  as  the  gradient  projection  technique.  The 
condition  to  be  imposed  on  the  controls  is  that  they  make 
the  system  identifiable,  that  is,  the  eigenvalues  of  M are 
nonzero.  Maximizing  the  weighted  trace  of  M is  not  mean- 
ingful for  this  case  since  there  will  always  be  sufficient 
information  in  the  output  to  estimate  I.  For  this  reason 
the  deterministic  case  is  not  addressed  further  in  this  re- 
search. 

If  measurement  noise  is  added  to  the  model,  then  Eqs 
(3)  and  (4)  become 

x(k*1)  = A(a)x(k)  * B(a)u(k)  (6) 

£(k)  = Cx(k)  * w(k)  (7) 

where  w(k)  is  an  r-dimensional  discrete  time  Gaussian  white 
noise  vector  with 

E(  w(k) ) = 0 and  E(w(k)w(j)T)  = R 6kj.  (8) 

Inserting  the  control  given  in  Eq  (1)  into  Eqs  (6)  and  (7) 
yields: 

x(k-*-1)  = [A(a)  + B(a)FyCjx(k)  ♦ B(a)v(k) 

♦ B(a)F  w(k)  (9) 
jr(k)  = Cx(k)  ♦ w(k)  (7) 

The  noise  term  w(k)  is  an  additional  uncertainty  beyond  that 
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of  a itself.  To  identify  a with  noise  in  the  measurements, 
it  is  desirable  to  have  the  output  £(k)  contain  as  much  in- 
formation about  a as  possible;  that  is,  the  weighted  trace 
of  M for  this  problem  should  be  as  large  as  possible.  It 
is  this  case,  measurement  noise  only,  that  is  addressed  in 
this  dissertation. 

If  measurement  and  process  noise  are  added  to  the  model 
then  Eqs  (3)  and  (Zf)  become: 


x(k+1 ) = A(a)x(k)  + B(a)u(k)  + *(k) 
£(k)  = Cx(k)  + w(k) 


(10) 

(ID 


Maximizing  a function  of  M is  meaningful  in  this  case, 
but  calculating  M becomes  very  difficult.  Mehra  (Ref  21) 
shows  that  calculating  the  conditional  likelihood  function 
In  p(  Yjj | a ) for  this  model  is  given  in  terms  Of  the  Kalman 
filter  model  corresponding  to  Eqs  (10)  and  (11).  The  Kalman 
filter  model,  as  presented  in  Mehra' s paper,  is 


&(kt-l)  = A(a)x(k)  «■  B(a)u(k)  ♦ K(k)y(k) 

— — o — — — o •—  — — 

v( k)  = Cjc(k)  ♦ Z(k)j/(k) 


where 


fc(k) 
*(k) 
£(k) 
P(k*1 ) 
z(k) 


E[x(k)|^(1),  2.(2) , • • • , ^(k-l )] 
[cP(k)V  ♦ g]**(i(k)  - Cfc(k)) 

A(a0)P(k)CT[cP(k)CT  ♦ r]”* 
A(a0)P(k)AT(aQ)  - K(k)S^(k)  ♦ £ 
[0P(k)CT  + gj* 
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(13) 
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The  measurement  and  process  noises  are  assumed  to  be  white 


and  Gaussian,  with 


E[f(k)]  = 0 

; E[w(k)j 

E[f(k)fT(  j)] 

= a*kj 

Efwtkjw1^  j)] 

= R«kj  • 

E[w(k)0T(  jj 

* 0 

EfxCOw^k)] 

= 0 

E[x(0)/(k)] 

= 0 

E[Kk)^T(  j} 

= i*kj 

= 0 


— » 


E[v(k)]  = 0 


As  shown  by  Mehra,  the  main  reason  for  the  complexity  of  M 
is  that  the  covariance  matrix  of  Z(k)v(k)  is  time  varying 
and  dependent  on  a in  a nonlinear  manner.  Taking  the  par- 
tials  of  this  matrix  with  respect  to  a becomes  difficult. 

Once  M is  calculated,  it  is  shown  in  the  next  chapter 
that  afj/aa  has  to  be  evaluated  for  the  gradient  algorithm. 
Taking  the  partial  of  M with  respect  to  a increases  the 
complexity  of  the  problem  significantly.  Solving  the  prob- 
lem of  determining  the  optimal  control  for  this  model  with 
process  noise  is  considered  beyond  the  scope  of  this  re- 
search and  is  not  addressed.  The  work  developed  in  this  re- 
search should  provide  insight  into  solving  this  more  diffi- 
cult, but  practically  significant,  problem. 


Objective  and  Organization 

The  objective  of  this  research  is  to  develop  an  algo- 
rithm that  will  calculate  the  constant  feedback  matrix  and 

open-loop  control  term  that  optimally  aids  in  identifying  a 
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set  of  unknown  parameters  in  a given  linear  system  model. 

The  constraints  on  the  system  are  an  energy  constraint  on 
the  open-loop  portion  of  the  controls  and  a restriction 
that  poles  of  the  feedback  system  must  remain  within  a pre- 
determined constraint  space.  This  space  depends  on  the  sta- 
bility and  responsiveness  required  of  the  system. 

To  achieve  this  objective,  an  algorithm  is  developed 
in  Chapter  II  that  considers  the  case  of  a system  with  a 
single  input  control  and  a single  unknown  parameter.  This 
reduces  the  complexity  of  the  problem  and  makes  it  easier  to 
understand  the  concepts  and  procedure  used  to  develop  the 
algorithm. 

To  help  in  understanding  the  solution,  a three  dimen- 
sional problem  is  selected,  and  the  optimal  control  is  cal- 
culated. This  is  performed  in  Chapter  III.  Although  there 
is  only  one  unknown  parameter  in  the  system,  its  influence 
is  dispersed  throughout  the  system  and  control  matrices  in 
a nonlinear  manner.  Cases  with  both  two  and  three  dimen- 
sional output  vectors  are  studied. 

An  estimator  is  also  designed  to  demonstrate  enhanced 
parameter  estimation  with  use  of  the  affine  control.  The 
calculated  optimal  affine  control  and  a noise  generator  are 
applied  to  the  example  case  to  generate  system  outputs.  The 
estimator  then  computes  the  estimate  of  the  unknc^rn  param- 
eter by  using  the  output  measurements  and  input  controls. 

The  results  of  the  estimator  for  different  noise  levels  and 


different  number  of  time  sequences  are  presented 


i 


I 


i 


i 


In  Chapter  IV  extensions  to  the  problem  addressed  in 
Chapter  II  are  presented.  For  these  extensions  it  can  be 
easily  seen  by  the  complexity  of  the  algorithms  that  the 
computational  burden  increases  substantially.  Chapter  V 
concludes  by  summarizing  the  contributions  made  in  this 
research  along  with  recommendations  for  future  research  in 
this  area. 


■I 

I 


II  Optimal  Feedback  Control  Design 


Introduction 

Consider  a discrete  system  assumed  to  be  modelled  ad- 
equately by 

x(k*1)  = A(I)x(k)  + B(a)u(k)  (14) 

£(k)  * Cx(k)  + w(k)  (15) 

where  x(k)  is  an  n-dimensional  state  vector,  u(k)  is  a sca- 
lar control,  £(k)  is  an  r-dimensional  output  vector,  and 
w(k)  is  an  r-dimensional  measurement  noise  vector.  A(a)  is 
an  nxn  time-invariant  plant  matrix  that  may  have  confluent 
eigenvalues  (i.e.  nondistinct ) , B(a)  is  an  nxl  time-invari- 
ant control  matrix,  and  C is  an  rxn  time-invariant  output 
matrix  that  is  independent  of  the  parameter  a.  The  effects 
of  unknown  scalar  parameter  a are  confined  to  A(a)  and  B(a) 
only,  and  aQ  is  a nominal  value  of  a.  Vector  w(k)  is  a 
Gaussian  white  noise  sequence  with 

E w(k)  =0  ; E w(k)wT(j)  = R <5kj  (16) 

The  initial  condition  x(0)  is  assumed  to  be  known  ex- 
actly. It  is  also  assumed  that  the  system  described  by  Eqs 
(14)  and  (15)  is  completely  controllable  and  complete]y  ob- 
servable for  all  a,  i.e., 

rank  Jj3(a),  A(a )B(a) , • • • , An“*  (a)B(a  )j  = n 

rank  [c_T,  AT(a)CT,  • • • , ( An-1  (a)  )TCTj  = n 


Without  this  assumption  there  may  be  modes  of  the  system 
that  cannot  be  controlled  to  the  desired  value  and  the  un- 
known parameter  may  not  be  observable  in  the  output. 

The  control  u(k)  is  restricted  to  having  a linear  out- 
put feedback  term  and  an  open-loop  external  control  term,  a 
member  of  the  affine  control  law.  Thus, 

u(k)  = Fyi(k)  «•  v(k)  (17) 

where  Fy  is  an  1 xr  time-invariant  feedback  matrix  and  v(k) 
is  a scalar  control.  An  energy  constraint  is  imposed  on  the 
external  control  term  v(k)  such  that 

N-1  2 

£ v (k)  * E (18) 

k=0 

where  N-1  is  the  last  instant  of  time  the  control  is  applied, 
i.e.  N = final  time  instant.  Without  a constraint  on  v(k) 
it  would  be  desirable  to  have  as  large  a signal  as  possible 
in  order  to  get  an  output  with  a large  signal  to  noise  ratio. 
Practically  this  would  not  be  possible  since  there  is  a lim- 
ited amount  of  energy  available  to  control  most  systems. 

The  restriction  on  the  feedback  matrix  Fy  is  that  the  eigen- 
values of  the  feedback  homogeneous  system 

x(k*1)  * [A(a)  ♦ B(a)FyC]x(k)  (19) 

be  members  of  a constraint  space  * that  is  predetermined  so 
as  to  ensure  desirable  stability  and  time  response  charac- 
teristics (see  page  53).  Therefore,  the  objective  of  the 
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research  is  to  determine  the  optimal  feedback  matrix  Fy  and 
optimal  external  control  sequence  vector  V,  where  this  vec- 
tor is  = jjv(Q),  v(1  ),•*•,  v(N-1  )J  such  that  the  best 
estimate  of  the  unknown  parameter  can  be  found  while  main- 
taining the  constraint  requirements. 

It  is  assumed  that  the  estimator  that  will  be  used  to 
identify  a will  be  in  the  class  of  efficient  or  asymptoti- 
cally efficient  estimators.  That  is,  they  are  asymptoti- 
cally unbiased,  minimum  variance  estimators  such  as  maxi- 
mum likelihood  estimators.  Using  this  assumption,  the  der- 
ivation of  the  optimal  Fy  and  V can  be  separated  from  the 
derivation  of  the  estimator  (Ref  23).  Maximization  of  the 
Fisher  information  matrix  is  the  criterion  used  to  calculate 
Fy  and  V.  For  the  case  of  the  single  unknown  parameter, 
this  matrix  becomes  a scalar.  For  any  unbiased  estimate 

A 

a of  the  unknown  parameter  a,  it  is  known  that  a lower  bound 
of  its  error  variance  is  given  by  the  Cramer-Rao  inequality 
(Ref  23) 

E [(a-l)2|  ij  * M”1  (20) 

where  a is  such  that 


a 


The  Fisher  information  scalar  is  given  by 


M 


E r ( 3 in  p(ym|5)\ 

yNL  l 3a  j 


2j 

J 


(21  ) 


(22) 
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Throughout  this  re- 


where  YN  = j^Ck):  k = 0,  1 , - * • , Nj 
port  all  partials  are  evaluated  at  a = aQ. 

For  an  asymptotically  efficient  estimator  the  lower 
bound  is  asymptotically  achieved;  therefore,  by  maximizing 
the  information  scalar,  the  estimation  error  is  decreased. 
The  F^r  and  V are  found  that  maximize  M. 

Information  Criterion 

Substituting  Eq  (17)  into  Eq  (14),  the  equations  for 
the  system  become: 

x(k+1)  = [A(a)  * B(a)FyCjx(k)  * 3(a)v(k) 

+ 3(a)F  w(k)  (23) 

•j 

^(k)  = Gx(k)  + w(k)  (21+) 

The  feedback  system  has  both  process  and  measurement  noise, 
but  they  are  the  same.  It  is  shown  that  this  simplifies  the 
means  of  calculating  the  information  scalar. 

In  order  to  solve  M,  p(Yjj|  a)  must  be  calculated.  Now, 
by  repeated  application  of  Bayes*  Rule, 

p(yn|5)  = p(i(n)|yn_1,S)p(ynJ  I) 

= p(l(N)|  Yn-1  ,a)p(y(N-l  )|  V2’S)p(YN-2l 
N 

= TT  p(y(  3)1  Y . . ,a)  (25) 

j=1  “ 

From  Eq  (24)  if  x(0)  and  ^(0)  are  known,  w(0)  can  be  com- 
puted. From  Eq  (23)  if  w(0),  x(0),  and  a are  known,  x ( 1 ) 
can  be  computed.  This  process  can  be  repeated  so  that  if 
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x(0),  a,  and  are  known,  then  = w(k):  k = 0,  I,'-*, 

j-lj  is  known  and  so  is  x(j).  With  this  information  the 


following  can  be  written: 


p(l(*i)|  Yj_!  ,a)  = p(^(  j)|  ,a) 

= p(^(  j)|  x(  j) ) 
E[x.(  j)|  x(  j)]  = Cx(  j ) 

E[(^(j)  - Cx(  j ) ) (jr(  j ) - Cx(j))TJ 


= E [w(  j )wT(  j )j 


Therefore , 


p(*<j>|t,  ,,!> . 7=4=  .^uJr'acj) 

J 1 J 2ir  R L 


where  n(j)  = £(j)  - Cx(j)  . 


p(yN|a)  = f[— Lr 

. .72wIH  L 


' N , 

In  p(YNja)  = Constant  - £ £ E n(j) 


(3D 


Therefore 


31n  p(YJa)  N 3nT(j)  , 
Jil = - V ■ R n(  j ) 

3 a ja 


3 In  p(YN|a) 


-T  / . 


3 a 


N N 3 n ( j ) , rn  , 3n(k) 

= Z Z ■H"'1n(j)r(k)r  — =r- 


J U\  k=1 


3 a 


3 a 


•T,  . ■ 


_ r 3 n ( j)  i m i 3 n ( j ) 

- Z ~7-_-  RH(3)I(j)£  — 

j = 1 3 a 3 a 


N-1  N 3nT(j)  . m 

+ 2Z  2 — S s.( J)5  (k)R’ 


3 = 1 k=j*l 


3 a 


3 2.(k) 
3 a 


(33) 


To  obtain  the  information  scalar,  the  expected  value  of  Eq 
(33)  is  computed.  The  expected  value  of  the  second  term  on 
the  right  side  of  Eq  (33)  is  zero.  This  can  be  seen  from 
the  following  equations: 


n(j)  = £(j)  - Cx( j ) = w(j) 
3 5 ( 3 ) g 

— = — = — rfeCj)  - Cx(j) 

3 a 3 a L J 


— C_ 


3 x(  j) 

3 a 


(34) 

(35) 


A terra  of  the  double  summation  is 


3hT(;j)  1.  i3n(k) 

2 = R n(j)nx(k)R  - — where  k > 3 

3a  — — — - g ~ 

T 

3 X ( j)  m 1 m 1 3 x(k) 

= 2— ~ — C IT  ‘ w(  j )£i  (k)if 1 C — — 

3a  ~ " 3a 

T 

m 1 Dx(k)  3 X ( j)  m 1 

= 2wA(k)H",C; = - — C1^  1(3) 

3 a 3 a 


(36) 


From  Eq  (23) 
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3 x(k)  3 x(k-1 ) 

= [A(a)  - 3(S)FyCj 


3a 


3 a 


* ~~~Z  r^Ca ) * B(a)F  c]x(k-l  ) 

3 a L ” 


3 B(a ) 33(a) 

-v(k-1)  - — — Fj(k.l)  (37) 


9 3 


as  -y- 


It  can  be  seen  that  3 x(k)/ 3 a is  a function  of  w(k-1), 
w(k-2),*--,  7(j),...,  w(0)  and  that  3x(  j)/3  a is  a func- 
tion of  v?(;j-l),  w( j— 2 ) , • • ♦ , w(0).  The  expected  value  of 
Eq  (36)  is 

3 x ( k ) 3XT(-) 


r .T  -.1  ° ± 

E " 


— v ' rp  1 -J 

I— c- lR  w(J)] 


9 a 


9 a 


r-T/,.0  -Tn-U  3-(k)  r>  - 


= 2£  [iT(k)J  e[r"  C- 


3 a 


— cV’wCj)]  (38) 


because  vf(k)  is  independent  of  the  other  terms.  Since 
E Fw( k )J  = 0,  Eq  (38)  is  zero. 

The  information  scalar  then  becomes: 


M = E 


N 3 nT(  j) 


Lj=i 


5 <3. 


.m  « 3 n(  j)  . 

■R  n(  j)nr(  j)R  1 I 

3 a J 


(39) 


The  summation  is  a scalar  and  using  the  fact  that  the  trace 
of  a scalar  equals  the  scalar  and 


tr  [XYZ]  = tr  [zxy]  = tr  [yzx] 
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where  X,  Y_,  and  Z,  are  matrices,  we  obtain,  after  substi- 
tuting Eqs  (3i+)  and  (35)  into  Eq  (39) 


" N r 3x(j)  3Xi(j)  rn  1 m 1 

M = E V tr  — ^ ~ j— O' lR“ '«(  j)wx  ( j)R  C 

j _ ^ 3 si  3 3 


Since  E 


N N 

Z trh  ■ 2 trE  [xj] 

Lj=1  J j=1 


and  it  was  shown  that 


3 x(j)/  3a  is  independent  of  w( j),  then 

N r 3 x(  j)  3 x ( j)  rn  , T 1 I 

M = ^ trE = C^R  w(j)wi(j)R  C; 

L 3 a 3a  J 

N I r 3X(  j ) 3 X ( j ) m 1 r<_  _rn  -1  _1 

= jrtrjEj 1 —CR  E[w(j)»l(j)j8  C (41) 

j=1  L L 3a  3 a J 

This  is  valid  because  [s  x(j)/  3a][3xT(j)/  3a]  is  independent 
of  w(j)wT(j)  and  the  expected  value  of  the  product  of  two 
independent  variables  is  equal  to  the  product  of  the  ex- 
pected values  of  the  variables.  Equation  (41 ) becomes 


N r r 3x(  j ) 3 x ( j )|  m 1 

M a y tr  E — — (TR  C 

j_1  3a  3 a J 


because  E [w(  j )wx  ( j j]  = R. 

Combining  Eqs  (23)  and  (24)  with  Eq  (37)  to  form  an 
augmented  system,  the  equations  become: 
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Using  this  notation,  Eq  (43)  becomes: 


XA(k+1)  = AAXA(k)  + BA(k)v(k)  + £Aw(k)  (48) 

c,  /x  f3  x( j)  3 xT( j) 

Since  an  expression  in  Eq  (42)  is  E , it  is 

3 a 3 a J 

necessary  to  compute  this  expected  value.  To  get  this  let 
E[xA(k*1)]  «XA(kH) 


(49) 


= AAXA(k)  + BAv(k) 


then 


3 x(k)  3 xT(k) 


3 a 


3 a 


= E 


hTXA(k)X^(k)h( 


= hT[XA(k)X^(k)  * PA(k)]  h (50) 


where 


“T  ■ [2  ! I]  n*2„ 

I is  an  nxn  identity  matrix 


(5D 


PAa=>  ■ AAPA(k-, )Aa  . DaHD^ 


(52) 


PA(0)  = 0 since  x(0)  and  3x(0)/  3 a = 0 and  are 
known  exactly. 


Now 


M = £ tr  hT  [XA(  j )xj(  j ) + PA(  j )]  hcV1  C 

3 = 1 L 


N 


'Pa  , Arp  T — 1 

= ^ tr|hlXA( j)X^( 'c 


3=1 


N 


. £tr 
3=1 


hTPA(  j )hCTR“1  c"j  (53) 


J 


The  term  inside  the  brackets  in  the  first  summation  is  a 
scalar.  Using  the  properties  of  a trace,  M becomes 


N 

z 

3=1 


N 


M = r j)  + Z^^V^hcV’c 


(54) 


3 = 1 
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If  there  were  no  feedback,  then  the  information  scalar 
as  given  by  Mehra  (Ref  20)  would  be: 


M 


Z j) 

j=1 


(55) 


The  differences  between  Eqs  (54)  and  (55)  are  that  in  Eq 
(54)  the  vector  XA(j)  is  a function  of  F^  while  in  Eq  (55) 
it  is  not,  and  that  Eq  (54)  includes  the  additional  term 
which  is  called  J(Fy).  This  additional  term 


J(V 


« m m i 

Z tr  & Ea^— - 

j=i  L 


(56) 


represents  the  weighted  trace  of  the  error  covariance  for 

the  estimate  of  the  vector  3 x(j)/  3a  where  the  weighting 
T -1 

matrix  is  C R C. 

From  Eq  (52)  it  can  be  easily  shown  that 

£a(1°  = 2 ^2 aSSa^a*1’1  (57) 

PA(0)  = 0 (56) 

From  Eq  (56)  it  can  be  seen  that  for  a large  uncertainty 
about  how  a change  in  the  parameter  a affects  a change  in 
the  state  x(j),  i.e.  large  ?A(k),  J(Fy)  is  increased.  We 
are  dealing  with  h PA(k)h,  the  lower  right  partition  of  the 
matrix  PA(k)  and  with  3x(k)/  3a.  This  means  that  more 
information  about  a is  gained  by  feeding  back  the  output  if 
J(Fy)  is  always  positive  which  is  true  and  proven  in  Appen- 
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dix  A.  If  the  uncertainty  is  small,  little  additional  in- 
formation is  gained  by  feedback.  From  Eqs  (56)  and  (57)  it 
can  be  seen  that  measurement  noise  increases  the  error  co- 
variance  matrix  but  decreases  the  weighting  matrix.  This 
implies  that  the  magnitude  of  the  measurement  noise  may  have 
a small  effect  on  this  additional  term.  For  example,  con- 
sider the  case  where  R = rl  . The  term  CTR_1C  becomes 
icTC  and  Eq  (57)  becomes: 

£a<«  * r 1 4*  (59) 

3=1 

Therefore  from  Eq  (56)  it  is  seen  that  the  r terms  cancel 
and  that  J(F  ) is  independent  of  the  measurement  noise. 

To  see  the  effect  of  the  feedback  vector  on  J(F,  ),  Eq 

y 

(57)  can  be  expanded  to 


1=1 


= F RF  T V A1-1B  B^A7)1-1 
-y— , y Z,  -A  —A— A ^ —A 

1=1 


(60) 


since  F RF  is  a scalar.  From  Eq  (56) 

«/  J 


N 


J(Fy)  = 2 t* 
>1 


h\RFyT  | A^B.B^A^^’hcV’c 

1=1 


= Z ^vTtr fh 1 c 


^ u-y— y - 
.1=1  L 


(61) 
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because 


N j 1-1  N N-i 

2 2 J*  - (62) 

j=1  1=1  j=1 

and 

tr^X(j)  = £ tr£(  0)  (63) 

j=l  j=1 

From  Eq  (61)  it  is  easily  seen  that  if  there  is  no  feedback, 

F = 0 , then  J(F  ) = 0.  Another  important  fact  about  J(Ft  ) 

is  that  this  term  can  never  be  negative  (proved  in  Appendix 
A).  This  implies  that  the  additional  term  can  never  de- 
crease information  about  the  parameter  a because  of  feed- 
back. 

The  expressions  for  the  information  scalar  and  the  con- 
straints are  known.  It  is  now  necessary  to  apply  an  opti- 
mization technique  to  maximize  the  information  available. 

This  is  done  by  using  a gradient  technique. 

Optimizing  the  Information 

The  optimal  F must  be  such  that  it  maximizes  the  sum 
of  the  two  terms  in  Eq  (54).  To  find  the  V and  F that  max- 
imize M,  a gradient  method  is  used.  The  reason  a first  or- 
der method  is  selected  is  that  the  complexity  increases 
significantly  if  higher  order  partial  derivatives  are  re- 
quired. This  will  become  more  obvious  subsequently  when  the 
first  order  partials  are  derived.  Another  advantage  is  that 
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the  gradient  method  will  always  converge  at  least  to  a lo- 
cal maximum,  whereas  some  higher  order  methods  only  converge 
if  the  initial  conditions  are  close  to  the  optimum  value. 
Disadvantages  are  that  the  global  maximum  may  not  be  ob- 
tained and  that  convergence  is  usually  slow  near  the  opti- 
mum value.  In  using  the  gradient  method,  the  gradient 


3M 


F „ 
-yo 

v 

-o 


has  to  be  evaluated.  This  would  require  taking  the  partial 
of  a scalar  with  respect  to  an  (r-^N)  dimensional  vector  (see 
Appendix  B for  an  explanation  of  this  procedure). 

Evaluation  of  this  partial  at  each  iteration  step  would 
consume  a large  amount  of  computation  time.  It  is  shown 
that  only  the  r-dimensional  vector 


3 M 


’V 


F 

-y 


v 


has  to  be  calculated.  In  order  to  do  this,  the  expression 
for  M can  be  expanded  so  as  to  simplify  the  partials. 

Modified  Form  of  M.  The  first  term  of  Eq  (54)  for  M 
can  be  written  again  as 

K = 2 ^(k)hCTR“1ChTXA(k)  (64) 

k=1 
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Proof  of  this  is  shown  in  Appendix  G. 
Also 


(63) 


.K~  - 

-A^A  " 


Ak  ( a ) B ( a ) 


(69) 


t-i-toa)B(a) 
j 3 a L F 


because  the  bottom  term  of  this  matrix  can  be  written  as 

3 3(a) 


_3 

3a  L' 


Ap(a)B(a) 


3 Ap(a)  k 0 ii 

-F_B(a)  + Aj(a) 


(70) 


3 a 


3 a 


and  it  is  the  right  hand  terms  that  are  obtained  when  AA  is 
postmultiplied  by  B^. 

Let 


(k)  = 
~aa 


g(0)f  v(0) , v(1  ),•  v(k-1)l 


(71 ) 


Elements  n*1  through  2n  of  Vx  (k)  equal  Zer°  becaU6e  °f  th® 
equation  3 x(0)/  3a  = 0 so  Eq  (65)  becomes 
Ak(a)  1 Ap”1  (a)B(a)  

1 3 B(S) 


xA00 


3 Ak(a) 


L 


3a 


i-rAk-1(a)B(I)]J 

3 3L'F  J| 


3 a 


l . 
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N 

K = 2)  V^OK)E(k)JaCTR""1CJhTET(lt)VxOE) 
k=1 

Let 

W(k)  = E ( k ) hC ^ R" 1 ChT E ^ ( k ) 

then 

N rn 

K = £ Vi(k)W(k)Vx(k) 
k=1 


(75) 


(76) 


(77) 


(78) 


It  is  shown  later  in  this  section  that  it  is  advantageous  to 


express  K in  a quadratic  form  without  the  summation.  This 
is  accomplished  by  expressing  K as 


K = VJ(N)WnVx(N) 


where 


v£(N)  = x T C 0 ) , v(0),  ▼(!),•••,  v(N-1 ) 


-(1 }(n+1 )x(n+1 


2(n-i  JxCn-ijJ  (n*N)x(n<.N) 


-(2)(n+2)x(n+2) 


— (N-2)x(!I-2)J  (n+N)x(n+N) 


* L^(N)j  (n+N)x(n*N)  (8l) 

Substituting  the  values  for  W(k)  into  Eq  (81)  yields  Eq  (82) 
which  is  shown  on  the  next  page.  If  x(0)  = 0 , the  first 

n columns  and  n rows  of  are  zero  and  can  be  deleted. 

r * 

Let  the  remaining  matrix  be  ^ and  the  elements  of  are 


Z — ^(ajBCa)  TCTR“1  C 


k=max(i, j) 


“ rte“j(a)B(S)] 


where  J(F  ) is  given  by  Eq  (61).  For  the  scalar  input  case, 
this  equation  can  be  further  simplified. 


tr 


jB  ab^  ( a*  )n“  -WV 1 c 


f~r 


= tr 


L 


, T.N-jg  ~|  [VT / .T \N-j.  ^Tq-1  r' 

h —A  IaJ  ^ S 1 


= ( A^ ) N“ J hCTR“ 1 ChT A^“ jBA 


(35) 


T T 

because  for  X,  Y both  n-dimensional  vectors,  trXY  = Y X ; 
therefore , 


N 


J(Fy)  “ Z RFyTB^ ( a\ ) N“ JhC V 1 ChTA:2“ 


■A 


(36) 


j=1 


To  maximize  M a straight  gradient  approach  could  be 
used,  but  as  mentioned  earlier,  this  would  be  very  diffi- 
cult. A simplified  method  can  be  used  because  of  the  special 
structure  of  Eq  (84)  and  is  discussed  in  the  next  subsection. 

Maximum  Value  of  M.  It  is  general  knowledge  that,  giv- 
T 

en  a quadratic  X AX  where  X is  an  n-dimensional  vector,  A 

T T 

is  an  nxn  matrix,  and  X X = 1 , the  maximum  value  of  X AX  is 

X , the  maximum  eigenvalue  of  A (Ref  22).  The  value  of  X 
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T 

that  gives  the  maximum  of  X AX  is  the  unit  eigenvector  cor- 
responding to  A . It  is  this  fact  that  is  used  to  simplify 

LU 

computing  the  maximum  of  M. 

From  the  work  leading  up  to  Eq  (84)  it  is  shown  that 
Wjj  and  J(F  ) are  functions  of  F only  and  not  of  the  input 
controls  while  V is  naturally  a function  of  only  the  input 
controls.  This  implies  that  for  a given  value  of  F , the 
matrix  and  scalar  J(F  ) are  determined.  Using  the  pre- 

1/ 

T 

sented  theory,  if  V V = 1 the  maximum  value  of  M would  be 

'W  = WV  * <87> 


where  Mm(F^)  is  the  maximum  value  of  M for  a given  F , 

X max(£y)  is  the  maximum  eigenvalue  of  for  a given  F . 

For  the  problem  considered  in  this  research,  the  con- 
straint on  the  controls  is  given  in  Eq  ( 1 8 ) which  is 


^ v2(k)  * E (18) 

k=0 


This  can  be  expressed  as 

VTV  * E (88) 

however,  in  order  to  use  the  above  approach  the  constraint 
on  the  external  controls  would  have  to  be 

VTV  * E (89) 

It  is  shown  in  Appendix  D that  for  the  cost  function  in- 
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volved,  Eq  (84),  all  the  energy  will  be  used  and  in  fact 

VTV  = E is  valid. 

T 

Since  V V equals  E instead  of  one’ 


T 

1 1 

v/F  s/i~ 


(90) 


W =Xmax(VE  + J(V  (93) 

The  value  of  the  control  sequence  V that  gives  this  maximum 

value  is  the  eigenvector  of  magnitude  JW of  corresponding 

to  the  eigenvalue  X „(F  ). 

max  — y 

The  above  work  shows  that  if  an  optimal  feedback  matrix 
had  been  determined,  finding  the  optimal  control  sequence  is 
relatively  easy.  The  difficult  step  is  to  determine  the 
optimal  F that  maximizes  M (F„)  while  constraining  the 
eigenvalues  or  poles  of  the  feedback  system  to  a given  con- 


straint space.  The  overall  optimal  information  scalar  will 


under  the  constraint  that  the  maximum  energy  is  E and  the 
eigenvalues  remain  in  4>. 

The  procedure  for  calculating  the  optimal  values  of  F 

J 

and  V can  be  restated  for  clarity.  Use  an  optimization  tech- 
nique, which  will  be  shown  in  the  next  section,  to  compute 

the  optimal  F„  and  it  will  be  called  F „ . . Once  F 4.  is 
-y  -yopt  -yopt 

1 

computed,  calculate  and  its  maximum  eigenvalue, 

^max^-yopt^*  Next  compute  the  eigenvector  of  correspond- 
ing to  Xmax(Fy0pt),  which  is  called  Vq  t,  and  scale  it  so 
T 

that  T£opt^opt  = ^ * The  vector  is  the  value 

for  the  control  sequence.  As  mentioned,  the  difficult  part 
is  finding  F ^ and  the  next  section  discusses  the  approach 
taken  to  obtain  this  Ixr  matrix. 

For  the  case  where  x(0)  / 0 , the  optimal  initial 
state  vector  is  computed  along  with  the  optimal  control  se- 
quence. A constraint  must  be  placed  on  the  initial  condi- 
tion vector  as  well  as  the  controls.  The  information  scalar 
for  this  case  is 


M * - J<£y> 


(95) 


and  assume  the  combined  constraint  on  is  ^ 


then 


W =Xn.ax<£y>E  *J<£y> 


(96) 


where  X „(F  ) is  now  the  maximum  eigenvalue  of  WM  and  Vv  is 
max  — y — w —A 

the  corresponding  eigenvector.  The  procedure  for  finding 
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M 


the  F to  maximize  X (F  ) that  will  be  developed  directly 
— y max  —y 

applies  to  this  situation.  The  optimal  Vx  then  contains 
both  the  external  control  sequence,  V,  and  initial  state 
vector,  x(0). 

It  is  also  possible  to  put  assigned  weighting  values  on 
the  elements  of  the  initial  state  vector  and  control  se- 
quence. Let 

E'  = hf*f(0)  * bfxfco)  h2,1v2(l)  * h2.2v2(2) 

*•••*  ‘>^n.,v2(»-D  (97) 

Then  if 


hn*N-1 


e'  = (|[Vx)T(5Vx) 


For  this  case 


M « (Evx)TWn(HVx)  ♦ J(F  ) 

* $r£VUx  * j<v 


(100) 


and  X max(£y)  is  now  the  maximum  eigenvalue  of  jfHjfjjff  and  Vx 
is  the  corresponding  eigenvalue. 


Gradient  Technique  For  Computing  F 

For  a given  F , the  maximum  value  of  M is 
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W * X m«(  VE  ‘ J<Ey> 


(93) 


A first  order  gradient  technique  is  used  to  find  the  F,  that 

S 

maximizes  Mffl(Fy).  This  requires  the  computation  of  the  par- 
tial 


3Mm<y 

— T 


T T 

F 1 = F x 
-y  -yc 


which  from  Eq  (93)  can  be  computed  as  follows: 


3M  (F  ) 
m -y 

* — 

8£y 


-y 


-yc 


3X  (F  ) 

= E 53S^=r- 

3J(Fv) 

* — =tf- 


3-y 


F 1 = F x 
-y  -yc 


F T = F T 
-y  -yc 


(101 ) 


where  F_  is  the  current  value  of  F . It  can  now  be  seen 
-yc  -y 

that,  as  mentioned  earlier,  the  partial  of  the  information 
scalar  with  respect  to  only  the  feedback  lxr  matrix  is  re- 
quired. This  significantly  reduces  the  computation  re- 
quired since  the  dimension  of  F is  much  smaller  than  the 

y 

dimension  of  V. 

s ( • ) 

The  partials  ■*  A are  evaluated  at  the  current  iter- 

•V 

ation  value  of  F so  the  notation  indicating  the  point  of 
evaluation  is  henceforth  eliminated.  The  partial  of 
Xmax(Fy)  with  respect  to  can  be  obtained  from 
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3X  N N ^ ma„(F)  3*,. 

max  — y ^ max  — y 1 


(102) 


where  w^  are  the  elements  of  Symmetry  can  be  used  to 

simplify  the  calculations  and  Eq  (102)  becomes: 


3X  max(V 


N 1-1  3X  (F  ) 3w. 

2Z  z ...max  -y ij 

i=1  j=l  3wij i *  3-y 

. I 3Xmax(V  3wi: 

Z 1 — J Tl 

l-l  swij 


(103) 


To  obtain 


31  _(F  ) 


, the  fact  that  is  real  and 


symmetric  is  used.  This  means  that  there  can  always  be 
found  an  orthogonal  transformation  matrix  T such  that 


is  * f*il 


(104) 


where  ^ is  a diagonal  NxN  matrix  with  components  being  the 

i 

eigenvalues,  x 's,  of  W^,  T is  an  NxN  orthogonal  matrix  such 
T -1 

that  T = T and  whose  columns  are  the  unit  eigenvectors, 

i 

e's,  of  Wjj.  Choosing  X1  to  be  \ max  gives 


^max 


= -- - j*u  1 • • • ' ®N 

• L l I 1 . 


(105) 


where  e^  is  the  normalized  eigenvector  corresponding  to  the 
eigenvalue  X^.  It  is  easily  seen  that 


^max^-y)  ~ ^max  -N^max 


(106 


Moreover,  Porter  and  Crossley  (Ref  26)  show  that 


3X  (F  ) 
max  — y ' 

T d~~ 


3 w.  . 
1J 


6 9 

max.  max. 


(107 


where  e 


max, 


is  the  k-th  element  of  e 


-max 


This  means  that 


the  eigenvector  corresponding  to  the  maximum  eigenvalue  of 


must  be  calculated  to  obtain  the  gradient. 


In  order  to  calculate  must  be  computed  first. 

-max  — N 

From  Eq  (83)  it  can  be  seen  that  is  symmetric  and  real. 
Also  from  Eq  (64)  it  can  be  seen  that  K ^ 0 , so  W,,  is 

positive  semidefinite. 


To  calculate  the 


n(n+1  ) 


oartials, 


3 

3a 


Ap"i(a)B(a) 


■ g ■ ■■  ■ independent  terms  in  the 
, have  to  be  evaluated  and  Aj,(a) 

is  given  by  Eq  (66).  This  computation  can  be  done  iter- 
atively. 


3 


where 


Ap(a)B(a)  = B(a) 

and 

3Ap(a)  3A(a)  3B(a) 

Z = 1 — + — 1 F C 

3a  3 a 3 a y“ 


Also 


(109) 

(110) 


' N 3 

w.  . = y — i. 

k=max(i, j) 

N 


A^-i(i)B(5) 


TCTR"  1 C— i — j ( a 


3a  L 


A£“ J(S)3(S) 


k=max(i+1 , j+1 r 


Ap~i(a)B(a) 

tgth~1  c — 

r J 

~ ~ ” 3 a 

\k- j /-sa.-ri 


* -^r[iF"i(i)2(a) 


W^c— 


A^“J(a)B(a)| 
3 alT 


k=max(i, j) 


= Vl,j+1 


3 a L 


A^“X(5)B(S) 


tctr“1g— 

3a 


'a^-^  (a)B(a  )j 


(111) 

k=max(i, j) 


Therefore,  first  calculate  the  elements  in  the  bottom  row  of 
and  then  use  Eq  (111)  to  obtain  the  other  elements  in  the 
lower  triangle  of  -N*  Knowing  the  elements  of  it  is  then 
possible  to  get  X max(Fy)  and  e^^,  the  unit  length  eigen- 
vector  of  by  means  of  existing  algorithms  referenced  in 
the  book  by  K.  Rektorys  (Ref  29). 


The 
Using  Eq 


f 


required  in  Eq  (102) 


(83) 


must  now  be  obtained. 
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and  from  Eq  (111)  the  remaining  terms  can  be  obtained  iter- 
atively by  means  of 


- ° 


013 


3F 

-y 


3 -y 


3 aL 


Ap”1  (a)B(a)]  VS’1  £ — 


3 -y 


~ A^(S)B(a) 


3 a L 


k=max(i, j) 


3 a L 


A^“j(S)B(5) 


Vs’1 2.—^ 


3 F. 


-y  l 


3 aL 


Ap”1 ( a ) B ( a ) 


k=max(i, j) 


where 


3F. 


“j4"m(5)s(S> 


-y 

derived.  From  Eq  (108) 


is  an  nxr  matrix  and  has  to  be 


O 


Pi  *“) 

— IAp(a)B(a) 

= Ap(a) -m 

lAp”1  (a)B(a)| 

_ 3a  L r “ _ 

J s£y 

_ 3 a L r J 

TifCa)—  fk?-1  (S)B(S>1 


SAp(a) 

3 a 


— :a5_1(S)B(S)|  (114) 

5 V L-  - i 


where 


aV 


3 

Ap  ( a ) 3 ( a )1 

_ 3a 

L“*  “ J 

3 F. 


-y 


3 B(a) 


3 a 


= 0 


(115) 


Using  the  proofs  in  Appendix  E,  the  following  equations  can 
readily  be  generated. 


■^-mATr(a)-^rA°"1(5)3(a) 

3Fy  F 3a  l F 


J 


= 3(a)- 


3a  L 


a?"1 (5)3(5) 

•"r  “ 


TCT 


(116) 


r 3 A-ri(a)~if 


3% 


iip 


3 a J 


Ap”1 (a)B(a) 


3 B(5) 


3 a 


Ap”1 (a)B(a) 


T T 


(117) 


3Ap(a)  * 


3 a 


3 F 

-y 


A?“1(a)B(a) 
; F ~ j 


a A.  (a)  ra-1  , , - T - 

F-  --  2 AJ  (5)S(S)S  (a) 


3 a 


1=1 


A*(a) 

— r 


m-1 -1CT 


(113) 


Equations  (116)  through  (118)  and  Eq  (108)  can  now  be  sub- 
stituted into  Eq  (114).  This  seems  to  be  complex,  but  the 
equations  are  iterative  and  can  be  calculated  rapidly  on  a 
computer. 


All  the  tools  are  now  available  to  calculate 


aWi 


since 


( N i-1  3X  (F  ) 3w 

2 r ' 

. , . , 3 w.  . 3 F 

3X  (F  ) 1=1  J=1  iJ  7 

max  -v  _ ) 


N 3*  -(Fv)  3 w . . 

+ V 

^ 3 w 3F  J’ 

i=t  wij  -y 


3X  ( F ) 

..  ma?  ~y  is  evaluated  via  Eq  (107)> 


0 03) 


3 w.  . 


H by  means  of  Eqs  (112)  and  ( 1 1 3 ) » 


— fa(a)B(a)|  by  Eqs  (108),  (H6)  through  (118) 

17b  s L-f  - JJ 


and  ( 1 1 k ) , 


— — A1?(a)B(a)  through  Eqs  (108)  to  (110). 

3 a L"F  “ - 

3J(F  ) • n 

The  next  step  is  to  calculate  ifr-.  Rewriting  Eq 


(86)  yields 


«».>  • i a-Vc^^Vst*^  <«> 


so 
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3 J(F  > N 


= z 2 V1  chV]- jsA 

j=1L 


= 2 2 JFyRjB^A^^^hcV1^^-3^ 


-A  ~A 


♦ F.  TbT(  a:  )N“  jhCTR“1  C— 2-* : hTA;'"jB 


lylJCA-)1'  JhC"R  'C—^r^A^  jBa 

3-y  - 


(119) 


From  Eq  (69)  and  the  definition  of  h,  Eq  (119)  becomes: 


J(EJ 


2 2 J£vR(  — fAp~j(a)B(a)lTCTR“1C 
_ , y"  35  _"F  " „ " " 


• ^rr4?“j(a)B(5)l  * F T^;rAS"j(5)3(5)lT 


*V  ••  —M 

^ 3 aL  * 


• cV’c — LT!-LrAp"j(5)B(5) 

3 £v  L3  al  F 


(120) 


Equations  (108), (109), (116)  through  (118),  and  (114)  are 

3 J(F  ) 

used  with  Eq  (120)  to  evaluate  ■—  4-  for  a given  F . It 

3 F 1 y 

-y 

can  easily  be  seen  from  Eq  (120)  that  with  no  feedback, 


3*1(0) 


= 0 . 


All  the  information  is  available  to  compute  a . 


The  following  equation 


021  ) 


i i 

H 


M = FT  + 
o 


a m (f  ) 

m 


3 F 

-y 


AFy  * H.O.T. 
F = F 

-y  rye 


'/there  H.O.T.  are  the  higher  order  terms  of  the  Taylor  expan- 
sion, AF  is  the  step  3ize  in  the  feedback  space,  and  FT  is 
the  cost  function  before  the  step  was  taken,  shows  that  the 
step  should  be  in  the  same  direction  as  the  gradient,  if 
there  are  no  constraints,  inorder  to  maximize  M.  There  are 
in  fact  constraints  so  the  standard  gradient  approach  to 
maximizing  M must  be  modified. 


Transfer  Gradient  into  Complex  Space.  The  constraints 
are  defined  in  the  complex  space  of  eigenvalues.  An  example 
of  a constraint  is 


< 


S 


for  0<  5<1  and  i=1  , 2,---,n 


(122) 


which  states  that  all  eigenvalues  must  be  within  a circle 
whose  maximum  radius  is  one . This  is  a requirement  for  the 
system  to  be  stable.  Adding  another  constraint  which  is 


|lmag  A I ^ real  A (123) 

0 ^ real  A < 1 (124) 


represents  a sector  as  shown  in  Figure  3*  The  shaded  region 
is  *.  This  constraint  represents  a system  which  must  remain 
stable  and  the  frequency  of  oscillation  of  the  response  is 
restricted  to  be  less  than  a desired  amount. 


J 


53 


j W 


Figure  3.  Constraint  Space  in  Z-Plane;  R = 1 

At  this  point  the  constraints  are  in  the  eigenvalue 
space  and  the  gradient  is  in  the  feedback  vector  space.  It 
would  be  difficult,  and  insight  would  be  lost,  if  one  tried 
to  transform  the  constraint  area  in  the  eigenvalue  space  to 
the  feedback  space.  Instead  it  is  proposed  to  transform 
the  gradient  into  the  complex  eigenvalue  space.  That  is, 

3 w 

find  where  X v is  the  n-dimensional  vector  of  eigen- 

3LF 

values  of  Ap(a)  and  Ap(a)  is  given  be  Eq  (66).  Now  the  pro- 
cedure is  to  move  in  the  direction  of  this  gradient  in  the 
complex  space  and  be  sure  the  eigenvalues  stay  within  the 
constraint  apace. 


5k 


■Vitn  this  approach  a X_p,  instead  of  F^,  will  be  found 
that  maximizes  M.  Once  _Xp  is  computed,  a transformation 
equation  is  needed  to  compute  the  F that  is  required  to 

J 

transform  the  open-loop  eigenvalues  of  A(a)  to  feedback 
eigenvalues  given  by  A_p. 

To  obtain  an  equation  relating  the  feedback  eigenvalues 
to  Fy,  an  approach  developed  by  Porter  (Ref  25)  will  be 
used.  A detailed  description  of  the  approach  is  presented 
in  Appendix  F.  This  method  allows  one  to  obtain  a closed 
form  equation  for  F_,  where  F^  is  a state  feedback  Ixn  ma- 
trix,  in  terms  of  _Xp.  State  feedback  means  that  the  output 
is  the  entire  state  vector.  For  the  single  input  case  this 
is  not  a significant  achievement,  however,  this  method  also 
applies  to  the  multi-input  case  where  the  state  feedback  is 

now  an  mxn  matrix.  Most  available  methods  that  solve  for  F 

— s 

for  the  multi-input  case  are  not  in  a closed  form  solution. 
Some  of  the  methods  as  given  by  Porter  and  Crossley  (Ref  26) 
are  minimum-gain  controllers,  prescribed-gain  controllers, 
and  multi-stage  controllers.  The  dyadic  form  controller 
would  give  a closed  form  solution,  but  the  feedback  matrix 
would  have  to  have  the  dyadic  form 

le  = (127) 


where  £ and  z are  respectively  mxl  and  nxl  vectors.  Forcing 
the  feedback  matrix  to  be  of  rank  one  would  restrict  the 
degrees  of  freedom  allowed  to  select  an  optimal  matrix. 


The  approach  that  is  being  used  provides  the  relation- 
ship between  a state  feedback  matrix,  Fg  and  the  output 
eigenvalues  and  not  the  relationship  between  an  output  feed- 
back matrix,  F and  the  eigenvalues.  Since  this  research  is 
«/ 


concerned 

with  F , another  equation  is 

J 

needed  to  obtain  F 

from  F . 

It  is  known  that 

— s 

Is  3 

£y£ 

(123) 

where  C is  the  rxn  output  matrix.  If  C were  an  nxn  matrix 
of  rank  n,  then  F would  be 

J 

ly  = ZSC~]  (129) 

However,  in  general  this  is  not  the  case,  and  for  this  re- 
search C”1  does  not  exist.  It  is  shown  later  that  the 
pseudoinverse  of  C is  used  and  that  additional  constraints 
of  the  eigenvalues  are  required  to  obtain  F from  F . 

— y — S 

Computing  from  will  be  developed  for  the  multi- 
input system.  Even  though  this  Chapter  deals  only  with  the 
single  input  case,  it  will  save  having  to  repeat  much  of  the 
procedure  for  the  multi-input  case  discussed  in  Chapter  IV. 
Wherever  a significant  change  occurs  for  the  single  input 
case,  it  will  be  mentioned  at  that  time. 

Computing  the  State  Feedback  Equation.  The  method  of 
computing  the  state  feedback  equation  is  to  transform  the 
original  state  equation 

x(k+1)  = A(a)x(k)  + B(a)u(k) 
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(130) 


I 


I 

l 


to  the  generalized  control  canonical  form 

z(k+1)  = F, (a)z(k)  ♦ G (a) w (131) 

■“G  "**  “"'C  *™c 

where  F (a)  is  block  diagonal  (F,  , F,  . . , F.  ) and  G (a) 

^ -k1  -k2  -Km  -c 

is  block  diagonal  ).  For  the  reader  who 

*12  m 

is  unfamiliar  with  the  generalized  control  canonical  form  it 
would  be  helpful  to  read  Appendix  F at  this  time.  The  Ap- 
pendix explains  the  procedure  described  here  and  also  gives 
the  reasoning  behind  each  step.  Eq  (130  is  obtained  from 
Eq  (130)  by  means  of  three  transformations: 

(1)  x(k)  = Tz(k)  where  T is  an  nxn  matrix  and  det  T 
does  not  equal  zero.  Equation  030)  becomes: 


z(k+1 ) = T‘1A(a)Tz(k)  + T”13(a)u(k)  (132) 

(2)  u(k)  = Zv  (k)  v/here  Z is  an  mxm  matrix  and  det  Z 

— — c — — 

does  not  equal  zero.  Equation  032)  becomes: 


z(k+ 1 ) = T“1A(a)Tz(k)  ♦ T"1B(a)Zv.(k) 


033) 


(3)  The  introduction  of  feedback  v (k)  = w (k) 

** 0 * "0 

- Lz(k)  where  w (k)  is  an  mxl  external  control  vector 
and  L is  an  mxn  feedback  matrix.  The  purpose  of  the  L 
matrix  is  to  add  the  necessary  freedom  to  change  the 
transformed  A(a)  matrix  to  the  generalized  control  ca- 
nonical form  F (a).  The  w (k)  control  vector  still 
— c — c 

allows  one  to  apply  a state  feedback  to  the  generalized 
form  and  that  will  be  demonstrated  later.  Equation 
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(133)  becomes: 


z(k+ 1 ) = 


(T“]A(a)T 


T”1 B(a)ZL)z(k) 


* T*"1 B(a)Zw  (k) 


and 


E^(a) 

2c(a) 


= F (a)z(k)  *■ 


a-1 


Gc(a)w,(k) 


= T"  (A(a)T  - 3(a)ZL) 
= 


(135) 

(13D 


(136) 

(137) 


Appendix  F explains  how  the  matrices  T,  Z,  and  L are 
derived  once  the  desired  form  has  been  determined.  In  order 
to  move  the  feedback  eigenvalues  to  the  desired  locations,  a 
feedback  matrix  must  be  calculated.  Let 


w (k)  = Hz(k) 

— c — 


(133) 


where  H is  the  mxn  feedback  matrix  that  feeds  the  general- 
ized state  vector  into  the  input.  This  will  move  the  eigen- 
values and  Eq  ( 1 3 1 ) becomes: 


z(k+1)  = (F  (a)  + G (a)H)z(k) 
= F^z (k) 


(139) 


where  F^  is  the  desired  feedback  system  matrix.  It  is  se- 

. » i 

lected  to  be  in  the  block  diagonal  form  (F^  , Fd  , • • ■ 

k1  k2 

• • 

, Fd  ) v.nere  _Fd  is  a matrix  block  which  has  the  same  di- 


al 


m 


mension  as  the  corresponding  block  in  F,  and  it  is  in  the 
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companion  form.  Matrix  H can  now  be  evaluated  and  as  shown 


in  Appendix  F,  its  form  is 


^11’  ^12*  » 


4 1 

I 2(^+1),  , h2(it1  -»-lc2)l 


(140) 


l^n^n-k^),'*',  hmn 


where  k1  ♦ k2  * * * ' + kffl  = n are  the  control  invariants  (see 
Appendix  F).  For  the  single  input  case 


H = h1 1 , h12>  • ‘ hiJ 


(14D 


The  nonzero  elements  in  the  i-th  row  of  Eq  (140)  and  the 
only  row  of  Eq  (141)  are  the  coefficients  of  the  character- 
istic polynomial  of  kj_  and  n feedback  system  eigenvalues 
respectively.  This  can  be  expressed  in  a summation  form  of 
the  products  of  these  eigenvalues. 


k.-l+l  k,-l*2 

k 1 — 1 V-1 


2 2 ' ' * 2 xf.xf.-  • * xf. 

v w1  H^i-r1  1 2 1 


(142) 


where  the  maximum  number  of  summations  is  equal  to  k and 

m 

the  number  of  eigenvalues  in  each  product  of  the  summation 
is  equal  to 


Z ki  “ p * 1 

i-=1 


043) 


Equations  (142)  and  (143)  show  that  as  one  computes  the  el- 
ements in  a given  row  and  moves  from  left  to  right,  the  i.,- 
tn  summation  is  removed  first  and  then  the  i^_^-th  summation. 
This  continues  until  only  the  i^  summation  remains.  For  the 
single  input  case,  j = 1 and  = n . An  example  will 
better  explain  this  procedure. 

Assume  k^  = 3 and  the  desired  feedback  system  eigen- 
values are  Xp  , Xp  , and  Xp  . The  characteristic  polynomial 
is 


0 


M *2 

From  Appendix  F it  is  shown  that  the  elements  in  H come  from 
the  bottom  rows  of  the  companion  matrices.  Therefore,  the 
elements  in  the  first  row  of  H for  this  example  are 


I 


hi  i = ^ f Xt?  X F 

h12  = -(XF,XP2  * XF,XFj  * XF2XFj> 

h13  * XF,  * XF2  *X  Fj 

hH  to  h,n  = 0 


Using  Eq  (142),  for  h.  1 


(145) 


(146) 


0 47) 


(148) 


1 . k,  - 1 ♦ 1 ■ 3 


hi  = (“1}  2 2 2 V Y V. 

il-1  W1  W1  1 2 3 

= X F1  XF2X  F5 


(149) 


for  h. 


1 = k.  - 2 ♦ 1 = 2 


h12  * (_1 ^ 2 2 X F.  X F. 

i1  =1  i2=i1  + 1 ^ 2 

3 V Y * Y XF  + Y X F ^ 
M J?2  if  3 *2  3 
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(150) 


and  for 


1 = k1  - 3 + 1 = 1 


*13  - ^ 2 V = XF,  +XF2  +XF, 

i1 =1  X1 


(151  ) 


Therefore,  it  was  shown  for  this  example  that  Eq  (142)  is 
valid.  The  matrix  H is  a function  of  the  feedback  system 


eigenvalues. 


In  order  to  return  to  the  original  state  space,  it  is 


observed  that 


u(k)  = Zvc(k) 


and  from  Eqs  (134)  and  (133) 
u(k)  = ZHz(k)  - ZLz(k) 

= ZHT**1  x(k)  - ZLT**1  x(k) 
= Z(H  - L)^”1 x(k) 


(152) 


Therefore,  the  feedback  matrix  is 


F = Z(H  - L)T" 

“C  * 


(153) 


This  is  the  desired  equation  and  a means  of  obtaining  Fy 

from  F is  now  addressed. 

— s 

Computing  the  Output  Feedback  Equation.  Repeating  Eq 
(128)  yields 


F = F C 

— s — y— 


(128) 


and  postmultiplying  by  C1  yields 


F CT  = F CCT 

-S-  -y — 


054) 


T T *1 

where  C£  is  of  rank  r.  Therefore,  (CC  ) exists  and 

Fy  = FsCT(CCT)_1  055) 


F T = (CCT)"]CF  T 0 56) 

J b 

T -1  + 

The  term  C (CC  ) is  called  the  pseudoinverse  of  C or  C . 

It  is  not  an  exact  inverse  in  that  even  though  Eq  0 55) 

gives  an  Fy  from  Fg  and  C,  using  this  F in  Eq  (126)  will 

generally  not  yield  the  original  F . In  this  case  CC+  = I, 

but  C C / I.  The  Fy  obtained  is  the  solution  that  minimizes 

the  norm  of  F_  - F C and  whose  own  norm  is  the  smallest  of 

— 0 —y— 

all  solutions. 

This  problem  occurs  because  G is  only  of  rank  r;  there 
are  more  equations  than  unknowns  in  Eq  (128).  This  implies 
that  there  may  not  be  any  F 's  that  will  satisfy  Eq  (128). 
Equation  0 55)  gives  an  approximate  solution  and  if  this 
solution  is  used,  the  eigenvalues  of  Ap(a)  may  not  actually 
be  and  may  even  fall  outside  of  the  constraint  space.  To 
correct  this,  the  scalar  equations  (128)  must  be  forced  to 
be  consistent.  Consistency  means  that  even  though  there  are 
more  equations  than  unknowns,  they  are  such  that  only  one 
unique  solution  will  satisfy  all  of  them.  7/ith  this  condi- 
tion of  consistency  imposed  on  Eq  (128),  the  F calculated 

J 


from  Eq  (155)  when  placed  in  Eq  (123)  would  yield  the  given 

F . As  will  be  seen,  forcing  "consistency"  will  inpose  more 

— s 

constraints  on  _Xp. 

For  the  single  input  case,  Eq  (123)  can  be  written  as 


I 

= ff  . f . • • f 

[V  V snJ 

LV  y2 

r 

1 1 * 

X 

c , 

c J . J c 

1 

—1 

-*2 1 | -n 

L 

1 1 J 

(157) 


Since  the  rank  of  C is  r,  there  are  r linear  independent 
columns  of  C.  7/ithout  loss  of  generality,  assume  the  first 

r columns  are  independent  and  the  matrix  consisting  of  only 

« 

these  columns  is  C . That  is 


C 


c 


i 

and  the  rank  of  C is  r.  Therefore, 


and 


(156) 


059) 


(160) 


There  may  be  an  optimal  set  of  independent  columns  that 
give  "better"  consistency  constraints.  How  to  select  this 
set,  if  it  in  fact  does  exist,  is  not  investigated  in  this 
research  and  is  suggested  for  future  research.  A suggestion 
is  to  select  those  columns  with  the  greatest  number  of  zeros 


in  it  in  order  to  simplify  the  computations.  The  remaining 

i 

n-r  columns  of  C are  linearly  dependent  on  the  columns  of  C 


c . = 

— x 


y k.  .c . 
^ ij-j 
d=i 


061) 


where  i = r+1  , r ’■2,*  • •«.  n . From  Eq  (157)  it  is  easily 
seen  that 


F c. 

-y-x 


= F yk.x. 
-y  Z,  ij-j 


= y k,  .F  c . 
Zj  Ij-y-.] 


" ^i/s. 

j=1  J 


062) 


The  k.  .'s  can  be  computed  from  C,  and  Eq  (162)  gives  the  re- 
x j ~ 

lationship  that  must  be  maintained  so  that  the  F computed 

*/ 

from  Eq  (155)  will  produce  the  original  F when  inserted  in 

■™G 

Eq  (128).  Since  the  elements  of  F are  in  terms  of  the 

"~S 

feedback  system  eigenvalues,  X_p,  as  shown  in  Eq  (153),  Eq 
(162)  is  then  added  constraints  on  the  feedback  system 


eigenvalues. 


Computing  the  Gradient.  To  find 


3M  (F  ) 


the  follow- 


65 


ing  equation  is  used 

(163) 

— F = -Fc 

where  denotes  complex  conjugate.  The  conjugate  is  re- 

quired because  the  inner  product  is  being  taken  in  the  com- 
plex eigenvalue  space.  For  the  single  input  case  that  is 

T 

3 F 

discussed  in  this  section,  — is  an  rxn  matrix.  The 

multi-input  control  case  will  be  discussed  in  Chapter  IV. 
From  Eq  (153)  and  the  fact  that  Z is  a scalar,  Z, 


8M  (F) 
m — y 

3 VP 

T* 

3F  1 

—y 

3^f 

3FyT 

\ - \ y 

-F  -Fc 

F = F „ 
-y  -yc 

3Lp 

3F  T , T3HT 

= (T-1)— — Z (164) 

3 *F  ~ 3 iF 


since  H is  the  only  matrix  in  Eq  (153)  that  depends  on  A_p. 
Therefore , 


3F  T T -1  -IT  3£T 

2L.  = (CC1)  ' C(T  'r Z 


(165) 


3iF 


3X_p 


3 H 


To  find 


■ , the  partial  of  each  element  of  H with 


respect  to  each  element  of  X.p  must  be  determined.  From 
Eq  (142) 
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k.-l+2  k,-l*3 

c-n1"1  ^z  JZ 


X_>  Jj' . 

il=1  i2=1l+1  1l-l=1l-2+1  1] 


•1  » i2’  ’ il-1 


X 7-1  • * -X  „ 

r . F . 


For  the  example  given  earlier 


(-1}  Z 2 


X _ X 


= X 


“ -(-1)  Z ^F. 

F2  it=1  X1 

V2 


i1 =1  i2=il+1  11  12 

ip  i-2  * 2 

3 


1 A3 


■(X™  + X™  ) 

F1  "3 


(166) 


12  = (-i)u  = l 

3Xf2 

3h, , 3 h- 

— lit  to  — la  = o 

3X-ri  3X-ri 

r2  F2 

which  is  easily  validated  by  taking  the  partials  of  Eqs 

(I45)  to  (148)  with  respect  to  X_  . Also 

*2 


T1 

3H" 

T 1 
3 H 1 | 

T 1 

3 HT 

3X_, 

3X  _ 1 

3X^  I 

3X_ 

— r 

TT  1 
*1  1 

T?  • 

2I  1 

F 

(167) 


Using  Eqs  (101),  (163),  (167),  and  the  conjugate  of  Eq  (165), 


can  be  obtained.  This  gives  the  direction  that  the 

32f 

feedback  system  eigenvalues  should  move  in  order  to  increase 
M.  If  equality  constraints  are  imposed  or  the  boundaries  of 
inequality  constraints  are  reached,  the  eigenvalues  must 
move  along  the  constraints  while  still  maximizing  M.  The 
gradient  and  value  of  M can  be  evaluated  at  each  step  and  a 
check  made  to  ensure  that  M is  indeed  increasing. 

To  handle  the  equality  and  inequality  constraints  there 
are  many  methods  that  can  be  used  (Ref  19).  One  method  is  tc 
use  the  equality  constraints  to  decrease  the  dimension  of 
the  problem.  This  is  a very  good  method,  but  in  general  is 
difficult  to  apply  if  the  constraint  equations  are  complex 
and  nonlinear.  Another  method  is  to  void  the  constraints  by 
means  of  penalty  functions.  This  method  adjoins  the  equal- 
ity constraints  to  the  cost  function  and  weights  each  con- 
straint. That  is 

Cost  = m - *T(x_,);=:*(x_)  ( 1 68 ) 

m — — — £ 

where  ^( _Xp)  is  the  vector  consisting  of  the  constraints,  and 
H is  a positive  definite  weighting  matrix,  with  iteratively 
increasing  eigenvalues.  Since  the  i;(X.p),s  are  supposed  to 
be  zero,  the  cost  function  penalizes  any  nonzero  value. 

This  approach,  hov/ever,  leads  to  a more  difficult  cost  max- 


imization  problem 


The  approach  that  is  used  is  the  gradient  projection 
method.  This  method  is  useful  for  both  equality  and  in- 
equality constraints.  The  basic  concept  is  that  a step  to 
change  ip  is  made  so  that  it  not  only  increases  the  cost 
the  most,  to  first  order,  but  reduces  the  equality  con- 
straints towards  zero  and  ignores  inequality  constraints  if 
ip  is  in  the  correct  region.  Information  on  the  gradient 
projection  method  can  be  found  in  Rosen's  articles  (Ref  30 
and  31 ) • 

The  inequality  constraints  that  define  * in  general 
would  be  simple  enough  to  handle  directly.  For  example  if 
* is  such  that  only  a certain  degree  of  stability  is  re- 
quired then  * could  be  chosen  as 


* S,  k » 


2,*  • *,  n;  E < 1 


(169) 


This  is  shown  in  Fig  if.  The  constraint  boundary  is  well 
defined , 


(Re  r * (Im  X-,  )‘ 


K2 


(170) 


and  it  is  not  difficult  to  keep  the  eigenvalues  on  this 
boundary  even  when  the  gradient  points  out  of  the  region. 

With  the  information  presented  in  the  previous  sub- 
sections, it  is  now  possible  to  outline  the  algorithm  that 
will  provide  the  optimal  F and  V that  maximizes  the  infor- 


Eigenvalue  Trajectories 


mation 


Algorithm  for  Solution 


The  algorithm  for  finding  the  optimal  feedback  matrix 
and  external  controls  is  shown  in  Figure  5.  The  algorithm 


is  initialized  with  the  values  F 


and  a„  as  shown  in 


block  1.  The  choices  of  Zy0>  2p0>  and  ao  are  naturally  not 
independent.  The  F must  be  such  that  the  eigenvalues  of 


An  initial  choice  of  F, 


X_p0  being  the  eigenvalues  of  A(a)  not  only  is  a reasonable 
choice,  but  also  guarantees  that  the  equality  constraints 
for  consistency,  Eqs  (128)  and  (162),  are  satisfied.  Any 


i 


Figure  5.  Algorithm  For  Obtaining  Fyopt  and  Kopt 

other  choice,  resulting  from  experience  for  example,  must 
also  satisfy  all  constraints  imposed  on  A.„. 


With  the  initial  values  determined,  block  2 shows  that 

I 

is  computed  from  Eq  (83).  Computing  X ( F ) and  e of 
1 max  — y -“'fnax 

t i 

as  outlined  in  block  3 can  be  performed  using  many  exist- 


ing methods.  The  one  used  in  this  research  was  EIGRF  sub- 


routine which  is  part  of  the  IMSL  pack  for  the  CDC-6600. 

Computing  J(F  ) from  Eq  (86)  and  then  using  Eq  (93),  the 

value  of  Mffi(Fy)  is  obtained.  This  is  shown  in  block  4. 

For  the  initial  pass  through  the  algorithm,  the  next 

step  would  be  decision  block  6.  For  other  cases,  decision 

block  5 is  to  check  to  see  if  M (F  ) had  increased  with  the 

m — y 

new  ip.  It  is  possible  that  the  incremental  step  taken, 

AXp,  may  be  too  large  and  cause  the  new  \v  to  decrease  the 
value  of  Mm(F  ).  If  this  does  occur,  the  procedure  (Ref  19) 
is  to  decrease  this  increment  step  and  repeat  the  Drevious 
steps.  This  is  shown  as  blocks  13  and  1i+. 

The  stopping  conditions  are  discussed  later;  for  now  it 


is  assumed  they  have  not  been  reached.  Since  the  new  has 

— T 

increased  Mm(F^)  it  is  now  necessary  to  select  another  step 
in  X_p.  The  partials  in  blocks  7 and  8 are  computed  from  Eqs 


(101),  (103),  (107),  (113),  and  (120). 
this  gradient  to  the  eigenvalue  space, 


In  order  to  transfer 


must  be  com- 


puted, as  referenced  in  block  9,  by  taking  the  complex  con- 


jugate of  Eq  (165).  Using  Eq  (163)  it  is  now  possible  to 
3M  (F  ) 

derive  2 — jL_.  This  is  shown  in  block  10. 

3 X-, 

— 1? 


73 


Knowing  the  gradient,  the  next  step  as  shown  in  block 
11  is  to  determine  the  increment  step,  ax_p,  that  will  be 
used  to  move  in  the  gradient  direction.  There  are  many 
ways  to  determine  the  size  of  AA.p>  but  the  method  found  most 
acceptable  in  this  research  is  to  use  a direct  control  on 
the  step  length.  This  means  that  jAX_p  | is  set  but  the  di- 
rection depends  on  the  gradient  and  constraints. 


With  AX-,  determined,  the  new  X-  becomes: 

r*  — 


new  = old  X-  + AX-^ 
— r -"i*  t 


(171  ) 


As  shown  in  block  12,  the  F^  that  produces  the  new  feedback 
system  eigenvalues  can  be  computed  from  Eqs  (153)  and  (155). 
The  complete  procedure  is  repeated  until  a stopping  condi- 
tion is  met. 

The  stopping  conditions  used  in  this  research  are 

(1)  When  the  change  in  M (F)  is  less  than  some  pre- 

m — y 

selected  value,  or 

(2)  When  a time  limit  has  been  reached. 

These  were  selected  to  limit  the  computation  time  while 

still  giving  good  results.  Ideally  one  would  want  to  stop 

only  when  the  maximum  value  of  M (F  ) is  reached:  this  is  not 

m — y 

practical  since  the  optimal  value  of  M (F  ) is  not  known. 

m — y 

Chapter  III  addresses  an  example  and  provides  results 
that  not  only  show  the  parameter  estimate  statistics,  but 
compares  this  to  the  estimate  statistics  in  the  case  of  ap- 
plying only  optimal  open-loop  controls. 


Ill  Optimal  Inputs  for  Three  Dimensional  Example 


Introduction 

In  this  chapter  the  developed  algorithm  is  applied  to 
an  example.  The  optimal  values  of  F^  and  V are  obtained  for 
many  different  conditions.  A maximum  likelihood  estimator 
is  then  designed  and  applied  to  the  example  to  estimate  the 
unknown  parameter  when  the  feedback  controls  are  used.  Re- 
sults from  the  experiment  are  tabulated  along  with  results 
when  using  an  optimal  open-loop  control  for  the  same  ex- 
ample. The  advantage  of  the  feedback  controls  is  especially 
obvious  wnen  the  level  of  the  measurement  noise  is  very 
large. 


Before  presenting  this  example,  the  maximum  likelihood 
estimator  is  developed.  The  estimator  algorithm  is  pre- 
sented along  with  equations  that  describe  the  estimator. 

Next  the  example  problem  is  given.  Using  this  example  system 
the  procedure  to  compute  and  VQpt  in  Chapter  II  is 

used  to  obtain  the  optimal  input  controls.  A nominal  value 
for  a is  used  in  obtaining  £y0pt  and  ^ However,  a dif- 
ferent value  of  a is  used  when  ^ and  V ^ are  applied  to 
the  system  to  get  the  output  measurements.  A noise  genera- 
tor is  used  to  corrupt  the  measurement. 

The  input  control  sequence  and  the  output  measurements 
are  then  inputs  to  the  estimator.  Results  of  the  estimates 
are  tabulated  and  compared.  Noise  levels  as  well  as  the 
length  of  the  control  sequence  V are  varied  for  evaluation. 
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A two  dimensional  output  system  as  well  as  a three  dimension- 
al output  system  is  studied. 


The  case  in  which  no  feedback  is  used  is  also  investi- 
gated. Using  the  optimal  open-loop  controls  only  provides 
a basis  for  comparison  with  the  feedback  system. 

Maximum  Likelihood  Estimator 

The  likelihood  function  for  this  problem  is  In  p(Y^ja). 
The  purpose  of  the  estimator  is  to  find  the  parameter  value 
a that  maximizes  the  likelihood  function.  This  is  equiva- 
lent to  finding  the  a that  maximizes  the  probability  that 
the  observed  measurement  values  resulted  from  the  unknown 
parameter  being  a.  Using  Eq  (31)  the  likelihood  function  is 


N 


L (a)  = Constant  - 1 £ nx ( j )R_1 n( j ) 


072) 


j = 1 


I _ 


Maximizing  L (a)  is  the  same  as  minimizing  L(a)  where 


N 


\T-~-l 


L(a)  = *£  “ £x(j)rR  (ZU)  " Cx(  j ) ) 

j=1 


(173) 


therefore , 


L(a)  = min  L(a) 
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To  solve  Eq  (174)  a combination  gradient/Newton- 
Raphson  method  is  used.  The  reason  that  the  combination 
method  is  used  is  that  the  Newton-Raphson  approach  converges 
faster  than  gradient  algorithms  near  the  solution,  but  may 
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example  is  the  direct  constraint  on  the  step  length,  since 


the  gradient  method  is  used  only  to  get  close  to  the  solu- 
tion. The  accuracy  and  complexity  of  the  local  parabola 
approach  is  not  needed  and  the  total  change  in  L(a)  is  un- 
known. Therefore,  k is  chosen  to  be 

’ s 


0 L(a) 


(177) 


The  .015  was  selected  from  experience  to  keep  the  number  of 

iterations  to  a reasonable  level.  For  this  problem,  the 

time  for  switching  to  Newton-Raphson  will  occur  at  the  first 

time  the  function  L(a)  starts  to  increase. 

From  Eqs  (175)  and  (176)  it  can  be  seen  that  -q-— — - 

3 a 

2j  / — v 

and  — must  be  evaluated.  From  Sq  (173) 


3 T ( g ) N 3 X(  j ) I 

~~  = - Z — S '(zU)  » CX(J)) 


(173) 


Z £ 


^ L' 
,i=i  LL 


2„/  . \ _T 


R (^( j)  - Cx( j)) 


r 3 x(  j ) ~]T  lf  3 X ( j ) — | 

C • — = R"  C ■ — 5 

L 3 a J L 3 a j 


3x( j ) 3 2x(j) 

To  calculate  and  * — 


, Eq  (14)  is  used. 


0 79) 


(180) 


3x(j)  3 A(a)  3x(j)  33(a) 

— - = z — x(j)  * A(a) ♦ u ( j ) 

3a  3a  3a  3a 


^xtj)  32A(a)  3A(a)  3x(j) 

— 12 — = Z2“x(j)  * 2 — I I 

3a  3a  3a  3a 


32x(j)  323(a) 

+ A(a) 15 — + * — u(j)  (181) 

3 a 3 a 


where 


^x(O)  3x(0) 


= x(0)  = 0 


(182) 


For  a given  problem,  A(a),  B(a), 


3A(a)  3 B(a)  3 2A(I) 

T > I > 12 

3a  3 a 3a 


3 2S(a) 


can  easily  be  computed.  With  this  information' 


and  a given  sequence  of  controls  V = v(j):  j=0,  1 , * - • , N— 1 , 

a given  a,  and  Eqs  (175)  through  (182),  the  iterative  steps 
for  a can  be  evaluated.  Figure  6 shows  in  block  form  the 

A 

algorithm  that  is  used  to  find  a,  the  estimate  of  a.  The 
next  step  is  to  select  an  example  and  evaluate  the  results. 


2 


3 


*■  L(a^) 

Less  Then  Pre- 
-Selected bl./ 

— ? — 


Figure  6.  Estimator  Algorithm 


20-1) 


♦ 0 u(  j ) (183) 


l(  j) 


1 0 0 


0 1 0 


x ( j ) ♦ w(  j ) 


for  the  two  dimensional  output  system  and 


*U> 


1 0 x(J)  ♦ i(j) 


084) 


085) 


for  the  three  dimensional  output  system.  This  particular 
form  was  chosen  for  two  reasons.  First,  the  system  is  not 
linear  in  the  parameter  and  second,  the  numbers  were  se- 
lected so  that  the  consistency  constraint  is  fairly  simple 
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in  order  to  keep  the  computational  burden  to  a minimum.  The 
two  and  three  dimensional  output  equations  were  chosen  so 
the  reader  sees  a case  where  C”1  exists  and  a case  where  C+ 
is  required. 

The  nominal  value  of  a is  a =1.0  . It  is  also  as- 

o 

sumed  that  the  initial  conditions  for  the  states  are  zero. 
The  measurement  noise  is  a Gaussian  white  sequence  with 


w( j)w(k)* 


(186) 


and  a Gaussian  noise  generator  is  used  in  the  example.  The 
covariance  matrix  R is  varied  for  different  runs  of  the  ex- 
periment. 

From  Eq  (183) 


(187) 


3 B(a) 
3 a 


32A(a) 

7 F” 


and 


-1 

0 


(188) 


[_-2e2(1-5j 

32B(a) 

n ■ are  then  derived  from  Eqs  (187)  and 
3 a 


(138) 
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Obtaining  Estimator  Inputs 

Inputs  that  are  required  by  the  estimator  are  the  control 
sequence  V and  measurements  made  on  the  system's  output. 
The  external  control  sequence  V will  be  and  will  be 

computed  using  the  algorithm  developed  in  Chapter  II.  The 
measurements  will  be  obtained  when  Vg  t and  Fy0pt  are  ap- 
plied to  the  system  with  the  true  a,  aT  = 1 .2  , while  the 

a used  in  the  filter  is  aQ  = 1.0  , and  the  output  is 

corrupted  by  a white  Gaussian  noise. 

For  all  cases  the  energy  constraint  on  the  external 
controls  is  unity.  The  length  of  the  control  sequence  is 
varied  from  two  to  eight  steps.  The  values  for  R are 
varied;  however,  the  R used  to  compute  F and  t is 
the  same  R used  to  generate  the  output  measurements  when 

— yopt  anci  lopt  are  aPPl^-ed  t0  t*16  system. 

There  are  constraints  on  the  eigenvalues  in  order  to 
give  the  desired  stability  and  time  response  character- 
istics. These  constraints  are 


< .9 


ReX,,  * 0 
Fk 


ImagXp  * Rex- 
*k  *k 


(189) 

(190) 

(191) 


These  inequality  constraints  determine  the  constraint  space 
For  the  two  dimensional  output  example  it  can  easily  be 
seen  that  the  rank  of  C is  less  than  the  number  of  states. 
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This  implies  that  an  equality  constraint  on  Xp  is  needed. 
Only  one  equality  constraint  is  needed  since  the  number  of 
states  is  three  and  the  rank  of  C is  two.  The  constraint 
will  be  computed  by  using  the  equations  developed  in  the 
previous  chapter  and  in  Appendix  F which  the  reader  should 
read  to  understand  the  following  example. 

With  a = 1.0  , the  system  equations  become: 


"o 

.13667 

0 

0 

x(j+1)  = 

-1 

0 

-.7733 

x(j)  * 

0 

0 

1 

1 »5_ 

1 

0 

1 


(192) 


(193) 


Since  the  controllability  matrix 


B(ao),  A(ao)B(ao),  A2(ao)3(aQ) 


S 


3 


0 

0 


1 


0 -.1057 
.7733  -1.160 
1.5  1.477 


(194) 


A 

which  is  rank  three  so  as  shown  in  Appendix  F,  £ = g. 
A«1 

Forming  £ gives 


7.317 

1.94 

1.0 

14.2 

-1.294 

0 

■9.465 

0 

0 

(195) 


so  L = [.205  -.91  1.5]  . Using  Eq  (142) 


H = 


A—  A„  A-r,  . .(An  At-i  * A ^ A *-1  ) » 

LF1F2F3  "l  "2  F1  F3  r2  a3 


(X 


F,  **,2*»,3)]  <20’> 


From  Eq  (153)  the  state  feedback  matrix  is 


Se  - 


(An  An  An  “ . 205  )>  (.91  — An,  Ani  “ Ap  A n, 

F1  F2  F3  F1  F2  F1  F3 


- Ap  An  ),  ( Ap  * Av  ♦ Av  - 1.5) 
F2  F3  F1  - 2 r3 


T"1  (202) 


Using  Eq  (162), 


fs3  = *31%  + k32fs2 


(203) 


and  from  C it  is  clear  that 


k31  = k32  = 0 


so 


fs  * 0 
S3 


(204) 

(205) 


With  this  constraint  on  the  eigenvalues  then 

Fy  *=  FsCT(C cV1  (155) 

will  give  the  desired  output  feedback  matrix.  Substituting 
the  value  of  C,  T~\  and  Eq  (202)  into  Eq  (155)  will  give  F 

y 

in  terms  of  Xp. 

The  algorithm  in  Figure  5 w as  used  and  part  of  the  re- 
sults are  tabulated  in  Tables  I and  II.  The  tables  give 


Table  I.  When  N = 8 


1 

Feed 

Iback 

No  Feedback 

Time 

R a .01 

R = 10. 

R = .01 

R = 10. 

0 

.625 

.640 

.378 

.378 

1 

.534 

.537 

.503 

.503 

2 

.425 

.418 

.529 

.529 

3 

.308 

.294 

.450 

.450 

4 

.194 

.181 

.307 

.307 

5 

.100 

.090 

.159 

.159 

6 

.036 

.031 

.051 

.051 

7 

-.001 

-.001 

-.002 

-.002 

only  the  optimal  open-loop,  Vq  portion  of  the  optimal 
feedback  control  along  with  the  optimal  open-loop  control 
if  there  were  no  feedback.  Table  I shows  both  controls  for 
the  case  when  N * 8 and  the  measurement  noise  levels  R = RI_ 
are  different.  Table  II  gives  the  open-loop  control  se- 
quences for  cases  when  the  total  number  of  time  steps  vary 
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from  two  to  seven.  It  can  be  seen  from  the  tables  that  with 


feedback,  most  of  the  energy  is  applied  earlier  in  the  se- 
quence. The  system  knows  that  the  feedback  controls  will 
provide  energy  later  in  the  sequence.  For  the  examples  used 
about  9^%  of  the  total  energy  came  from  the  feedback  controls 
while  about  6%  came  from  the  open-loop  controls. 

The  optima1  feedback  matrix,  FyQpt,  for  all  cases  ex- 
cept when  both  the  control  sequence  consisted  of  eight  con- 
trols and  the  diagonal  elements  of  R were  .01  is 


F 

— yopt 


1.940,  -.449 


(208) 


For  the  case  with  the  eight  step  control  sequence 


F 

-yopt 


jj  .173, 


362 


(209) 


Figures  7 and  8 show  the  trajectory  of  the  eigenvalues  while 
is  being  optimized.  Double  poles  occur  at  .75  and  .3  re- 
spectively because  of  the  constraint  that  the  sum  of  the 
eigenvalues  must  be  equal  to  1.5. 

For  the  three  dimensional  output  example,  no  equality 
constraint  is  required  since  the  rank  of  C is  three.  Table 
III  lists  VQpt  for  the  feedback  and  open-loop  cases  where 
the  control  sequence  is  set  at  six  steps.  The  optimal  feed- 
back matrix  is 


£-3.406,  1.966,  1.200 
The  eigenvalue  trajectory  is  shown  in  Figure  9. 


F 

-yopt 


(210) 
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Table  III.  ^or  Three  Dimensional  Output 


Time 

Feedback  No  Feedback 

N = 6 

N = 6 

0 

.756 

.553 

1 

.564 

.588 

2 

.309 

.483 

3 

.120 

.303 

4 

.023 

.131 

5 

-.004 

.019 

Figure  9.  Trajectories  for  Three  Dimensional  Output  Case 

The  feedback  matrices  and  optimal  external  controls  are 
used  along  with  a Gaussian  generator  to  produce  measurements 
at  each  time  step.  The  equations  used  are 


x(j*1)  * A(aT)  * B(aT)FyoptC  x( j) 

+ B(aT)v(j)  ♦ B(aT)Fyoptw( j)  (211) 
£(  j)  = Cx( j)  ♦ w(  j)  (212) 

where  a^  is  the  true  value  of  a which  is  1 .2.  The  actual 
measurements  along  with  the  controls  are  used  in  the  algo- 
rithm of  Figure  6 to  estimate  the  unknown  parameter  and  the 
next  section  discusses  the  results. 

Estimate  Statistics 

A Monte  Carlo  analysis  is  used,  composed  of  50  simula- 
tion runs  with  the  noise  generator  being  the  only  variable 
changing  on  each  run.  A statistical  mean  and  variance  of 
the  estimate  is  then  computed.  Tables  IV,  V,  and  VI  tabu- 
late the  estimates  for  the  optimal  feedback  and  optimal  open- 
loop  control  cases.  Tables  IV  and  V are  for  the  two  dimen- 
sional output  example  and  Table  VI  is  for  the  three  dimen- 
sional output  example. 


Table  IV.  Average  Estimate  and  Variance  For 


Table  VI.  Average  Estimate  and  Variance  For 
Three  Dimensional  Output,  N = 6 


Noise 

Feedback 

Open-loop 

R 

A A | 

avga  = a 

• — 0 

avg(aT-a)^ 

A | 

avga  = a 

• — 2 

avg(aT-a  r 

.10 

1 .2004 

.00019 

1.1988 

.00166 

1.0 

1.1985 

.00030 

1.2403 

.03013 

10. 

1.2016 

.00021 

1.2395 

.10035 

1.2018 

.00023 

1 .2833 

.12116 

Table  IV  shows  that  with  feedback  the  estimate  of  the 
parameter  does  not  vary  much  for  different  levels  of  noise, 
but  the  confidence  is  higher  for  the  cases  with  low  noise. 
For  the  open-loop  case,  both  the  parameter  estimate  and 
confidence  level  worsen  as  the  noise  level  increases.  This 
is  an  important  point;  there  is  a significant  performance 
benefit  in  using  feedback  over  open-loop  controls  when  R is 
large. 

Table  V shows  that  as  the  number  of  control  steps  in- 
creases, the  parameter  estimates  improve  and  so  do  the  con- 
fidence levels  in  these  estimates  for  the  feedback  control 
example.  For  the  open-loop  case,  this  improvement  is  not 
obvious.  As  displayed  in  the  table,  the  increased  number  of 
steps  should  improve  the  estimate,  since  more  data  is  avail- 
able for  the  estimator  and  more  inputs  are  allowed  to  help 
the  estimator.  As  should  be  expected,  the  table  also  shows 
that  the  estimates  are  better  for  the  lower  noise  levels. 
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Table  VI  basically  displays  the  same  results  as  Table  IV 
except  that  it  can  be  seen  that  the  confidence  levels  are 
better  mainly  due  to  the  faot  that  more  output  data  is  fed 
back  to  the  input. 

Figures  10  and  11  present  the  information  derived  from 
Tables  IV  and  VI  in  graphical  form.  It  is  seen  that  the 


.01  .10  1.0  10.  100 


.01  .10  1.0  10.  100 

Noise  Level  - R 


Figure  1 1 . Average  Estimate  ± One  Standard  Deviation 
Vs  Noise  Level  For  The  Three  Dimensional  Outnut  Eramnla 
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feedback  controls  improve  the  estimates  considerably  when 
there  is  much  measurement  noise  present.  For  low  noise 
levels  the  optimal  feedback  and  open-loop  controls  produce 
about  the  same  results.  If  there  is  no  noise,  both  controls 
should  produce  a nearly  perfect  estimate.  Figure  10  shows 
the  average  estimate  with  a ± one  sigma  deviation  about  it 
for  the  two  dimensional  output  example  and  Figure  1 1 is  the 
same  for  the  three  dimensional  example. 

Besides  the  fact  that  the  more  accurate  estimates  and 
variances  are  in  general  better  for  the  feedback  controls, 
another  performance  comparison  is  the  computer  execution 
time  required  by  the  estimator.  All  computation  was  done 
on  a CDC  6600  computer.  Table  VII  tabulates  this  information 
and  clearly  shows  the  time  saved  for  cases  with  high  noise 
levels  when  feedback  is  used.  It  is  understood  that  just 
giving  the  computation  times  is  not  conclusive  that  one  sys- 
tem is  superior  to  another.  Additions  and  multiplications 
required  per  method  would  be  better  as  well  as  the  number  of 
iterations  required.  However,  total  time  at  least  gives  an 
idea  of  the  relative  complexity  of  the  two  approaches. 

Conclusion 

The  results  of  the  previous  tables  and  figures  show 
that  for  most  cases,  the  feedback  controls  aided  the  esti- 
mator much  more  than  the  optimal  open-loop  controls.  As 
was  mentioned,  if  there  is  no  noise  present,  then  both  sys- 
tems should  give  the  same  results.  As  the  noise  level  in- 
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Table  VII.  Computer  Execution  Time 
Required  by  Estimator 


Steps 

N 


2.922 
2.891 
5.01 1 
4.429 
7.063 
6.494 


creases,  the  feedback  method  shows  a significant  improve- 
ment. For  small  noise  values  it  can  be  seen  that  in  some 
cases  the  open-loop  controls  actually  resulted  in  a better 
estimate;  however,  the  feedback  case  was  very  close  and  the 
accuracy  lost  by  using  feedback  control  was  small. 


IV  Multi-parameter  and  Multi-input  Cases 


Introduction 

In  Chapter  II  the  problem  of  using  feedback  controls  to 
aid  the  identification  of  a linear  system  model  was  solved 
for  the  case  of  a single  unknown  parameter  and  a single  in- 
put control.  The  purpose  was  to  simplify  the  problem  so 
that  a better  understanding  of  the  approach  to  the  solution 
and  development  of  the  algorithm  could  be  gained.  It  is 
shown  in  this  chapter  that,  for  the  more  complex  problems, 
the  basic  approach  taken  for  the  single  parameter,  single 
input  system  can  be  used. 

The  next  section  will  address  the  more  realistic  prob- 
lem of  how  to  compute  the  feedback  gains  and  external  open- 
loop  controls  when  there  is  more  than  one  unknown  parame- 
ter. The  following  section  looks  at  the  case  in  which 
multi-inputs  are  required.  Both  cases  increase  the  compu- 
tational burden,  and  it  may  be  practical  to  compute  input 
controls  only  for  small  dimensional  systems. 


Multiple  Unknown  Parameters 

With  multiple  unknown  parameters,  the  information  ma- 
trix of  Eq  (22)  becomes: 


In  p(YN|  a)~]r 

3 In  p(YN|af] 

Da 

3 a 

(213) 


where  a is  a p-dimensional  vector  of  unknown  parameters. 
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From  Eq  (32) 


For  the  same  reason  as  outlined  in  Chapter  II,  Eq  (£1 5) 
becomes : 


and  it  can  be  seen  from  Eqs  (34)  to  (42)  that  the  elements 
of  M become: 


100 


N r f3x(j)  3 XX(  j)l  T i 
mik  • I tr  E — “ ' =T—  c'R  'c 

it,  L L »*t  3akJ 


(217) 


£ (k*1 ) * xT(k*1 ), 


3 xT(k+ 1 ) 


Si  * ST(a) , 


3 BT  ( a ) 


T • 

rnrn  m 3 B (s) 

T7  1 ~ 


St  = F B^a),  F — 

i S7  7 3 ai 


A(a)  * B(a)FyC 


-A,  * 3A(a)  3 B(a)  I 

1 ;=-  * — F C A(a)  ♦ B(a)F  C 


taL  3 aL 


X.  (k+1 ) = A.  X.  (k)  * B.  v(k) 
”Ai  11  1 


Then 


W Am  rp  1 fTlA 

«u*  Z2l<J>££2  aV« 

J.1  1 


(218) 


(219) 


(220) 


(221  ) 


where  h1  is  defined  by  Eq  (51)  and 

PA  (J)  * A SA  (J-DaT  + D SSl 

Aik  A1  Alk  Ak  A1  Ak 

P.  (0)  - 0 
Alk 


(222) 


V tr[~hTP.  ( j)hCTR“1cl  (223) 

L“  “Aik  J 


(224) 


(225) 


Part  of  Eq  (222)  Is  not  affected  by  different  l's  and  need 
not  be  recomputed  for  each  i. 

The  dimension  of  the  information  matrix  is  pxp  so  a 
scalar  measure  of  M is  required  for  the  optimization  crite- 
rion. Chapter  I discussed  some  of  the  different  scalar 
measures  that  could  be  used.  These  criteria  include: 

( 1 ) maximize  tr  M 

(2)  minimize  tr  M”1 

(3)  minimize  jM"1) 

(4)  maximize  tr  (SM)  for  S a given  weighting  matrix 

It  would  be  desirable  to  use  the  second  or  third  cri- 

i terion  since  they  would  minimize  the  sum  of  the  parameter 

error  variances  or  the  volume  in  equiprobabllity  contours 
respectively.  Unfortunately,  obtaining  and  working  with  M”1 
is  very  difficult  and  mathematically  impractical.  For  these 
reasons  criteria  2 and  3 are  usually  not  used. 

The  first  criterion  is  the  easiest  to  work  with,  and 
for  this  reason  it  is  used  in  many  cases.  The  problem  with 
using  this  criterion,  as  stated  by  Zarrop  and  Goodwin  (Ref 
36),  is  that  it  may  lead  to  an  almost  singular  information 


matrix  which  has  large  diagonal  terms.  This  would  mean  that 


the  dispersion  matrix,  M”1 , would  also  have  large  diagonal 
terms,  that  is,  the  lower  bound  on  the  error  variances  of 
the  parameters  would  be  large.  This  problem  does  not  always 
occur,  but  is  something  that  must  be  carefully  checked. 

The  last  criterion  can  be  used  in  two  ways.  One  method 
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is  to  preselect  the  weighting  matrix  based  on  known  infor- 
mation and  then  hold  it  constant  throughout  the  algorithm. 
The  same  problem  that  can  occur  with  tr  M could  occur  here 
with  tr  SM,  so  it  is  not  recommended.  The  second  method  is 


to  select  S such  that  by  maximizing  tr  SM,  you  are  actually 


minimizing  tr  M”1  or  M-1  . An  algorithm  to  select  thin  S 


properly  is  presented  by  Gupta  and  Hall  (Ref  11).  Unfortu- 
nately the  derivation  of  S to  this  date  requires  that  the 
elements  of  M,  m^,  only  be  affected  by  the  input  controls. 
This  is  not  the  case  here  since  the  feedback  matrix  affects 
and  this  matrix  is  being  modified  in  the  process  of 
finding  £y0pt*  Attempts  of  extending  the  derivation  to  this 
case  were  not  fruitful,  so  this  eliminates  using  only  tr  SM 
as  a criterion.  It  is  easy  to  see  that  each  criterion  has 
its  advantages  and  disadvantages  and  a criterion  that  is  a 
combination  may  be  the  best.  For  this  reason  the  following 
ad  hoc  method  is  selected. 


The  criterion  to  maximize  tr  M is  used  with  the  addi- 


tion  that  the  value  of  M will  be  computed  at  each  itera- 


tion step.  The  overall  objective  is  to  decrease  the  volume 


in  equiprobablllty  contours,  i.e.,  decrease  M 


-II 


so 


.-1 


M | is  also  computed  at  each  step.  This  means  that  as  F 


-II 


is  being  iterated  to  maximize  tr  M,  M | is  being  monitored 


to  be  sure  it  is  decreasing.  If  a point  is  reached  where 


M | starts  to  increase,  F is  then  held  constant  and  the 

J 


criterion  to  maximize  tr  SM  is  used  where  S is  selected  such 
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trix  is  computed  using  the  algorithm  of  Figure  5,  while  the 
new  optimal  control  sequence,  Y.new>  is  computed  using  the 
method  outlined  in  the  weighted  trace  criterion  section. 

As  mentioned,  this  is  an  ad  hoc  method  and  certainly 
does  not  guarantee  that  the  global  optimum  values  for  F 
and  V are  found.  For  the  nonlinear  M derived  in  this  re- 
search, none  of  the  individual  criteria  given  in  this  sec- 
tion can  guarantee  a global  optimum  value  for  Fy  and  V. 

This  method  decreases  |m"^|  while  eliminating  the  computa- 
tional burden  required  when  working  directly  with  |m~^| . 

In  most  cases  the  criterion  of  maximizing  tr  SM  would  prob- 
ably not  be  required,  which  means  the  most  tractable  cri- 
terion would  be  used  throughout  the  algorithm.  This  occur- 
red in  the  examples  that  were  tested,  and  the  results  appear 
later  in  this  chapter.  The  next  subsection  will  develop  the 
max  tr  M and  max  tr  3M  criteria. 

Trace  M Criterion.  The  trace  of  M is  written  as 


tr  M » 


P f T ' ' / 

T.  rw  v ♦ j4(f.) 
L i 


• J \ J. 


(226) 


where 


(227) 
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w 


N 

2 


3 — T i BCa)"1  TGTR“1  ■ 


j=max(k, 1)  3aiL 


4s- 


i u 


3ai 


(228) 


N 


J-(Fy)  3 2 )N“JhCTR“1ChTA^“jB 

j=l 


i l 


x 'Ai 


(229) 


and  N is  the  number  of  external  control  steps  and  p is  the 
dimension  of  the  unknown  parameter  vector.  Equation  (226) 
becomes : 


tr  M = VJ 


W ’ + W ' + 

L-N1  -Ha 


I 1 

■ • + JLt  v 

j]%>  * 4%)  * 


rp  ti  ii 

• rusi  - J %) 


where 


ii  P 


3 2 — N. 

i=1  x 


3 2 J1(V 

i=1 


+ J (F  ) 
p -y 


(230) 


(231  ) 


(232) 


Equation  (230)  is  structured  the  same  as  Eq  (84)  so  the  pro- 
cedure used  to  maximize  M in  Chapter  II  can  be  used  here  to 
maximize  tr  M.  That  is,  for  F = F, 


■y  -yopt 


max  tr  M = max  VTWMV  * j"(F  ) 

V " “N”  W 
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(2 33) 


= X E + J (F  ) 


max 


-y 


where  X 


is  the  maximum  eigenvalue  of  W,t.  The  F 


max 


-N* 


M 

selected  so  that  the  maximum  eigenvalue  of  W„  is  X , and 

— N p ’ 

it 

V is  then  the  eigenvector  of  corresponding  to  X 


The 

V 

-yopt 

is 

is  X 

•q  } 

and 

Fmax 

to  X 

pmax 

and 

scaled  because  of  the  energy  constraint  E.  The  computation 

it 

has  increased  because  of  the  complexity  of  W^.  The  next 
subsection  develops  the  criterion  of  maximizing  tr  SM.  As 
mentioned  earlier,  this  criterion  is  required  if 
to  increase  while  tr  M is  being  decreased. 


M”1  starts 


Weighted  Trace  Criterion.  The  algorithm  discussed  in 
this  subsection  to  compute  V.0Dt  to  maximize  tr  SM  is  pre- 
sented by  Gupta  and  Hall  (Ref  11).  In  their  work  feedback 
is  not  considered,  but  if  feedback  is  used  and  held  constant 


their  approach  can  be  applied.  The  value  assigned  to  Fy 

I -1 1 

would  be  the  last  value  of  F before  |M 


started  to  in- 

-y  ■-  ■ 

crease . 

The  algorithm  steps  are: 

(1)  Let  be  the  value  of  V just  before 
creased. 


M”1  in- 


(2)  Let  M be  the  value  of  M before  M 

— o — 


-II 


increased. 


(3)  Find  a V with  energy  E that  maximizes  tr 

— m ^ — *o  — 


— m 


-o 


T -1 

and  that  yields  V V ^ 0.  M is  an  initial  weighting 

— m ~o 


matrix  and  M is  a function  of  only  V since  F has  been 

y 

determined.  If  tr  (M^MCV  ))  equals  tr  (M_1M(V  )), 

— — — m — — — n 
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stop,  since  is  the  optimal  value  for  V. 


(4)  Let  V,  = PV  + e v so  the  information  matrix  for 

— k — o — m 


-k  13 

5<k  - <■  2/3cMon 


(234) 


where  M is  the  information  matrix  for  input  V and  M 

— m — m — om 

is  the  "cross"  information  matrix  for  inputs  V and  V 

— o — m 

and  is  given  as 


rp  tl  It 

M = vHv„V  * J (F  ) 
— om  — o— N— m — y 


(235) 


The  energy  constraint  for  is  £ E which  re- 
quires that  and  e2VTV  «■  20eVTV  ^ E or  that 

—o—o  — m— m — o— m 


— o- m 


02E  ♦ €2E  + 2/3evTv  ^ E 

— o— m 


(236) 


thereby  constraining  the  choice  of  0 and  e. 

(5)  Now  use  Eqs  (234)  and  (236)  to  find  an  e between 


0 and  1 which  optimizes 


(c  is  chosen  arbitrarily 


to  start,  then  iterated  until  the  optimum  is  found; 
since  (234)  is  quadratic,  such  an  e will  eventually  be 
found).  This  step  requires  the  largest  amount  of  com- 
puter time  for  this  algorithm.  Once  e is  found,  0 can 


be  found  from  Eq  (236).  Gupta  and  Hall  show  that  if  V 


— o 


is  not  optimal,  then  an  improvement  can  always  be  made 

by  using  the  computed  0 and  e in  V,  = /3V  «■  t V . 

— k — o — m 

(6)  Check  for  termination  of  the  algorithm.  Termina- 
tion can  occur  when 
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a)  the  information  matrix  does  not  change  sub- 
stantially or 

b)  the  value  of  e is  very  close  to  zero. 

If  termination  does  not  occur  then  becomes  V^, 
becomes  M and  the  process  is  repeated. 

If  the  criterion  is  to  minimize  tr  M”1  then  in  step  3 
a V is  found  that  maximizes  tr  ( MT 1 MM" 1 ) and  in  step  5 an 
e is  found  that  minimizes  tr  . The  rest  of  the  procedure 
remains  the  same.  This  completes  the  criteria  derivation 
and  it  would  be  advantageous  to  present  an  example  at  this 
time. 


Two  Dimensional  Unknown  Parameter  Example.  The  example 
problem  is  designed  so  that  it  is  similar  to  the  example  in 
Chapter  III.  This  should  help  the  reader  in  following  the 
problem.  The  estimator  developed  in  Chapter  III  is  modified 
slightly,  as  will  be  shown  in  this  subsection,  in  order  to 
handle  the  two  dimensional  case. 

The  equations  that  model  the  example  system  are 


x( ) 


0 

-2  a, 

C 

1-a, 


.1367S| 


e 1 -2a £ 


-.1 


x(  j) 


109 


The  structure  of  A(a)  and  3(a)  is  chosen  so  that  when  the 
value  of  a^  is  set  to  jj.O,  .5J  , A(aQ)  and  BCa^)  are  the 
same  as  in  Chapter  III. 

The  initial  criterion  that  is  used  is  to  maximize  tr  M. 
At  the  same  time  M-1]  is  monitored.  If  at  anytime  during 
the  iterations  to  maximize  tr  M,  jM””1  j increases,  the  cri- 
terion to  maximize  tr  SM  is  selected. 

The  constraint  space  on  the  eigenvalues  is  the  same  as 
in  Chapter  III,  that  is 


x p < • 9 

*k 

(189) 

Rei,  ^ 0 

Tk 

(190) 

IraagXp  ^ RexF 

(191) 

The  values  of  the  covariance  matrix  for  the  noise,  R,  is 
chosen  to  be  R = 1.01  . 

The  algorithm  of  Figure  5 is  used  with  the  exception 
that  because  of  the  added  dimension  of  a,  the  partial  deriv- 
atives of  a variable  with  respect  to  a becomes  twice  as  com- 
plex. For  example,  with  one  unknown  parameter 
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E < 


I i 


! 


3»'  3 3T ( a ) 

■f  * 


3F 

-y 


3 a. 


CTR“1C- — — 


1 

Q-> 

L 35i 

L-'  " Jj 

li  -y 

and  with  two  unknown  parameters  it  becomes: 


(239) 


3w 


3V 


iLi  . 


3 a 


1C  3 

3 

Ap“^(a)B(a) 

- . „ T 

•s  TT 

3Fy 

c*  a^ 

_ 

3 ST(a)  m 

r*  -L  t-v  ' 


R C- 


3£y-L3a2 


-2-^"j(S)B(S) 


(240) 


During  the  iteration  steps  tr  M increased  from  577.4  to 
12,806.6  and  |m“^  decreased  from  .00001358  to  .00000016. 
At  no  time  did  |m“^|  increase,  so  the  weighted  trace  cri- 
terion is  never  needed.  The  optimal  values  for  F and  V 

V 


are 


F 

-yopt 

VT  : 
-opt 


1.948,  -.449 
.710,  .557,  .375,  .200,  .069,  -.002 


(241  ) 
(242) 


F . is  the  same  as  for  the  two  dimensional  output  example 
-yopt 

in  Chapter  III  when  R = 1.01  . This  is  not  surprising 

since  A(a  ) and  B(a  ) have  the  same  value  for  both  examples. 

— — o — — o 

The  partial  derivatives  given  in  Chapter  III  for  the  estima- 
tor are  computed.  Not  only  are  partials  required  as  shown 

32A(a)  32A(a) 

in  Chapter  III,  but  now  — and  — - — are  required  as 


3 a- 


l13  a2 


h 


in 


well.  These  are  needed  because  the  matrix 


i 


/ 


3 2L(|) 

3 2L(a) 

3 a2 

3 Sj3a2 

3 2L(a) 

32L(a) 

3 

3 a2 

2L(a) 

replaces  n in  the  Newrton-Raphson  algorithm. 
3 a 

value  of  the  parameters  is  selected  to  be 


The  true 


(243) 


The  modified  estimator  is  used  to  perform  the  Monte  Carlo 
analysis.  Estimates  are  made  with  the  feedback,  controls  and 
with  only  open-loop  controls,  as  in  Chapter  III.  The  number 
of  steps  in  the  control  sequence  is  set  at  6.  Table  VIII 
compares  the  results. 

Table  VIII.  Average  Estimate  and  Variance 
For  Parameter  Vector 


Feedback 

Open-loop 

avgl1  = a1 

1.19512 

1.22307 

A ^ f 

avgl2  = a2 

.40605 

.38544 

A p 

avg(a1T  - a1 ) 

.001 140 

.016953 

A 2 

avg(a2T  - a2) 

.000239 

.000450 

- 


i 
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The  optimal  open-loop  control  sequence  for  this  example  is 


VT  = f.621,  .612,  .437,  .214,  .054,  -.004 


(244) 


It  can  be  seen  from  Table  VIII  that  both  the  estimate  and 
the  confidence  in  the  estimate  are  improved  for  the  feed- 
back case. 

Even  though  the  weighted  trace  algorithm  to  minimize 
Iff1!  is  not  required,  it  is  used  to  find  the  optimal  open- 
loop  controls  to  minimize  |m“^|  . After  the  optimal  feed- 
back matrix  is  computed,  the  weighted  trace  algorithm  is 
applied,  using  this  feedback,  to  determine  if  a better  op- 
timal external  control  sequence  could  be  found.  Table  IX 
shows  the  results  using  both  the  max  £tr  mJ  and  max  |tr  SMj 
criteria  to  compute  the  open-loop  controls  and  Table  X 
shows  the  results  when  the  optimal  feedback  matrix  is  used. 


Table  X.  For  Feedback  Example 


Criterion 

Value 

tr  M 

tr  SM 

tr  M 

12806.6 

12779.5 

Ir'l 

.000000164 

.000000163 

t = 1 

.710 

.690 

2 

.557 

.511 

3 

.375 

.419 

4 

.200 

.275 

5 

.069 

.107 

6 

-.002 

-.005 

It  can  be  seen  that  when  the  weighted  trace  criterion 
is  used  |m“^|  improved,  but  very  slightly.  The  percentage 
improvement  is  well  below  1%  for  both  cases.  It  is  obvi- 
ous for  this  example  that  the  decreased  value  of  ) M*”1  ] ob- 
tained when  using  tr  SM  does  not  justify  the  additional 
computation  required  to  compute  S. 

In  summary,  max  tr  M criterion  is  used  because  the 
weighted  trace  criterion  can  not  handle  the  problem  in  which 
the  elements  of  M are  dependent  on  a variable  other  than  the 
controls  V.  The  problem  with  the  unweighted  trace  is  that 
M may  approach  singularity  and  cause  erroneous  results, 
thereby  necessitating  the  addition  of  a procedure  to  check 
for  this  condition.  In  the  algorithm  of  Figure  12,  this 
check  is  the  computation  of  Im”1)  . 


If  this  condition  occurs 


f t 

!l 


furtner  refinement  is  obtained  by  holding  Fv  constant  and 
varying  V only. 

The  modifications  to  the  basic  algorithm  of  Figure  5 to 
allow  for  multiple  parameters  are  listed  below.  The  other 


Table  XI.  Changes  to  Blocks  of  Algorithm  in 


Figure  5 to  Compute  VQ?t  and  £y0pt  ^or 
Multiple  Parameters 


Algorithm 
Block  Number 


Eq  Number  For 

Corresponding  Pesult 

Single  Unknown 

For  Multiple  Unknown 

Parameter 

Parameters 

33 

228,  231 

86,  87 

37,  229,  232 

- 

Algorithm  Fig.  12 

93 

Algorithm  Fig.  12 

blocks  are  only  affected  by  the  fact  that  the  computation 
is  increased  because  of  the  added  unknown  parameters.  An 
example  of  this  is  shown  in  Eqs  (239)  and  (24 0). 


Multiple  Input  Controls 

Chapter  II  addresses  the  single  input  control  problem. 
There  are  many  practical  applications  in  which  the  dimen- 
sion of  the  controls  is  more  than  one  and  this  section  mod- 
ifies the  previous  results  to  handle  this  situation. 

The  system  equations  (14)  and  (15)  become: 

x(k+1)  = A(a)x(k)  + 3(a)u(k) 

£(k)  = Cx(k)  * w(k) 
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(245) 

(246) 


gm 


1 ' 

where  3(a)  is  an  nxn  time-invariant  control  matrix  and  u(k) 
is  now  an  m-dimensional  control  vector.  The  feedback  control 
is 

u(k)  = Fy^rOc)  ♦ v(k)  (247) 

where  F is  an  mxr  time-invariant  feedback  matrix  and  v(k) 

~y  ~ 

is  an  m-dimensional  external  control  vector. 

A relationship  between  the  eigenvalues  of  the  feedback 
system  and  the  feedback  matrix  is  required.  This  relation- 
ship is  developed  in  Chapter  II  and  by  combining  Eqs  (153) 
and  (155)  is 

Fy  = Z(H  - L)T”1CT(CCT)“1  (£43) 

The  derivation  of  T,  L,  and  Z is  developed  in  Appendix  F and 
Eqs  (140),  (142),  and  (143)  show  how  H is  obtained.  The 
consistency  constraints  are  developed  for  the  single  input 
case  in  Chapter  II,  Eq  (162).  For  the  multiple  input  prob- 
lem, the  complexity  of  the  constraint  equations  increases. 

Let 


le  * £y2  (249) 

where  Fg  is  mxn,  Fy  is  mxr,  and  C is  rxn.  This  can  be  writ- 
ten as 


Again  since  C is  of  rank  r,  there  are  r independent  columns 
and  it  is  assumed  that  they  are  the  first  r columns.  Then 


f 

f • • • f 

T"  f , • • • f w 

S1 1 

S1 2 s1r 

y11  y12  y1r 

f 

f 

S21 

• 

= 

y21 

f ’ 

• * * f„ 

f‘  • • • f, 

_sm! 

emr 

ym1  ymr 

where 


- * [-1  j —2  I’  * ' [ — r] 

Repeating  Eq  (161) 


, * -1 


(251  ) 


(25 2) 


i * 


-i  = 2 kiJ^J 


(161) 


where  i = r+1 , r*2,---,  n and  k..  are  coefficients.  From 
Eq  (250) 


f = F c . 

-s-i  *y-i 


(253) 


where  f is  the  i-th  column  of  F . Combining  Eqs  (161) 

-s.i  s 


and  (253) 


1 


I 

I 

I 


r 


= yk..f 

1J-6. 

j'1 


3 


(254) 


Changing  this  vector  equation  to  m scalar  equations  gives 


f = y k.  .f 

“S1i  £ 1J"su 


f_  = y k.  .f 

^ t3“"sp  i 
3 = 1 


— s 


2i 


(255) 

(256) 


Is 


mi 


r 

y k,  .f 

<4.  i J— s . 
j = 1 mj 


(257) 


It  can  be  seen  that  since  the  dimension  of  the  controls  is 
now  m,  there  are  m times  as  many  consistency  constraint 
equality  equations  on  the  feedback  system  eigenvalues. 

Having  the  relationship  between  and  Fy,  Eq  (248),  and  the 
consistency  constraints,  Eqs  (255)  through  (257),  the  changes 
to  the  algorithm  in  Figure  5 because  of  multiple  inputs  are 
developed. 


Algorithm  Modification.  The  criterion  is  the  maximiz- 
ation of  the  information  scalar  because  it  is  assumed  that 
there  is  only  one  unknown  parameter.  If  there  were  multiple 
parameters,  the  material  developed  in  the  previous  section 
could  be  used.  The  information  scalar  M is  given  by 
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M=  Z SjCJ^V^Ch^C  j)  + Z tr  j )hCTR*”1  C (54) 

j = 1 j=1  L ’ "J 

where  the  relations  corresponding  to  Eq  (44)  to  Eq  (60)  are 


!A(j*1)  * ^A(j)  + 4^(k) 


ba  = BT(a), 


3 BT(a) 


3 a jmx2n 


(258) 


(259) 


r ~ A(a)  + B(a)F  C 

- - -y- 


-A  = 3 A(a)  3B(a) 


_ F C A(a)  ♦ B ( a ) F C 
D a J I J 2nx2n 


(260) 


7~|  TP  DTP  In!  r a T \ 1—1 


1=1 


(261  ) 


Equation  (54)  now  becomes: 


M = 2 fA(  j )hCTR“1  ChTXA(  j ) 


+ Z jtr  hTAA"jBAFyRFyTB^(AA)N“jhCTR"1C  (262) 
j=l  L J 


The  trace  operation  cannot  be  eliminated  as  is  done  in  the 
single  input  case.  It  is  proven  in  Appendix  A that  the  sec- 
ond term  is  non-negative. 

In  Chapter  II  the  procedure  was  to  express  Eq  (54)  into 


the  format 
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Equation  (84)  can  then  be  maximized  by  using  the  given  pro- 
cedures. The  following  equations  follow  the  work  given  by 
Eqs  (54)  through  (84),  but  are  modified  because  of  the  m- 

dimensional  control  vector.  The  increased  dimension  of  '!Lr 

—A 

and  V clearly  indicate  the  additional  computation  required 
to  solve  the  multi-input  problem. 

Equation  (65)  becomes,  for  the  multi-input  case. 


XAU)  = A£Xa(0)  + A^Bav(0)  - A^“2Bav(1) 


+ • • • ♦ AaB^v( j-1 ) 


—A’  -A  -A*  * ‘ -A 


2nx(2n+km) 


xA(°) 

v(0) 


(263) 


, K j-i ), 

— 1 (2n+km)x1 


Letting 


v£(j)  - 


xT(0),  vT ( 0 ) , vT(1),- vT(j-1) 


1x(n*km)(264) 


then  Eq  (73)  becomes: 


ET(k)  = 


4(i) 


Ap”1 (a)B(a) 


— -1 
r ip"1  (a)fi(aj| 


3 Ap(a)  | ^ 


3 a 


I B(a) 

.) 

I 3 3(a) 


3 a 


(265) 

2nx(n+km) 
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Eqs  (75)  and  (76)  remain  the  same  and  the  analogy  to  Eq  (78) 


K = Z v£oowooyx(k) 

k=1 


= V <V  V 


(266) 


where 


$ - 


xT(0),  vT(0),  vT(1),---,  IT(N-1)l1x(n,Nm)  (267) 


,Vhen  the  initial  condition  for  the  state  is  zero,  Eq  (266) 
becomes : 


T ' 

K = VxO 

— — N— 


(268) 


where  VT  = vT(0),  vT(1  ),•••,  vT(N-1 ) -lxNm  and  the  dimen- 

L i 

sion  of  V and  has  been  increased  by  a factor  of  m,  the 
dimension  of  the  input  controls.  This  significantly  in- 
creases the  computation  time  required  to  compute  (now 
a square  matrix  of  dimension  Nm  instead  of  N)  and  the  max- 
imum eigenvalue  of  and  its  corresponding  eigenvector. 
Equation  (82)  can  still  be  used  to  compute  since  the  only 
change  would  be  that  B(a)  is  now  a matrix  instead  of  a vec- 
tor. Using  Eq  (83)  to  compute  the  elements  of  can  be 
done;  however,  the  result  is  no  longer  a scalar,  but  as 
shown,  is  an  mxm  matrix  relationship. 

w^l  = £ -2-  Ap“k(a)B(a)  VtT’c  . 

i=max(k,  1)  3aL  -1 

’ 121 


(269) 


• -~Ap"’1(a)B(a) 


Following  the  procedure  of  Chapter  II,  once  the  elements 
of  are  computed,  the  maximum  eigenvalue  and  corresponding 
eigenvector  can  be  evaluated.  The  second  term  of  Eq  (84)  is 


J(F  ) 


z jtr 

j=i  L 


(270) 


3X  (F  ) 3w,  , 

To  compute  the  gradient,  the  ~ f and  ■ ■ ^ must 

9 2k!  3 V 

3X  (F  ) 

be  computed.  Computing  ^ is  not  difficult  since 

3 1 


Eq  (107)  gives  the  partial  of x raax(Fy ) with  respect  to  each 

I 

element  of  and  by  block  partitioning  the  result  becomes 


jAmaxv-v J 

r—r" " 

3 -itl* 


. Finding  the  other  partial  is  not  as  easy  since 


it  is  taking  the  partial  of  a matrix  with  respect  to  another 
matrix.  The  procedure  used  in  this  research  is  to  compute 

f 

-kl  for  each  value  of  k,l,a,  and  b where 


b = 1 , 2, • • •,  r 

k,  1,  a = 1,  2, • 


and  f.  is  an  element  of  F . Using  this  procedure 
ab  *’*•' 


i=max(k.,l)L  ° yab 


■ ^[Ap“k(a)B(a)]  T cV’c 

_ 3 3 - 


. -^[A^CaJBta)]  * -2=[A|-k(S)B(D]  Crg-'C 


3f  L3  a 
yab 


-~rAp"1(a)B(a> 


(271) 


where  the  iterative  equation 


_J L.  a^CDbCS)  = Av(a) 2 

3a  UF  - J - 3t 

^ -1  Ja.b 


— fA™-1  (a)B(S)l 
3 a[_  F J 


.[ap(S)|  — fe-1  (a)B(a)l  + — (~--”AF(a) 

L"F  J aaL*  J 3f  Laa  J 


. rif’CSJBCS)]  + 


a Ap ( a ) ^ 

3 a 3f„ 


."A““1(a)B(a)l  (272) 


-i-[A°(a)B(a)  = 0 

3 srF  J 


(273) 


can  be  used.  Equations  (271)  through  (273)  are  used  to  com- 

i 

pute  ~^1.-  and  by  performing  this  for  each  element  of  Fy, 
yab 


axmax^Fy^ 


can  be  computed. 
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It  would  be  useful  at  this  time  to  consider  an  example 
3X  (F  ) 

to  demonstrate  how  ■ ■ max.p is  computed.  Let  m = 2,  k = 

3F  1 


1 S 1,  SO 


Sima  x<£,> 

t ■nf—  = 

3 Ell 


-y 


3Xmax(V 
3 ( W1 1 ) 1 1 

5A~max^-v') 
3 (wn  )21 


3Xmax^-y ^ 
3(*11  )l2 


^max^-y^ 

3(wn  )22l 


(274) 


If  a,  b = 1 

in 

Eq  (271 ) then 

\ 

3 (wn  )n 

3 (*„  )12 

• 

’ill 

3 fY 
y11 

3fy11 

w 

3 ^1  1 

3 (wn  )2i 

3 (wn  )22 

* 

3fyn 

* P 

° y11 

(275) 


The  overall  result  that  is  desired  is  to  find  the  value  of 
3X  max^v5 

— — y i which  is 


3 f 


y11 


3X  max^-y^  _ J,  |, 

3 fy  L ^ 

y1 1 


3Xmax(-y)  3(v71 1 


i-1  J=1  3(wH)ij  3fyn 

* summation  for  (k,l)  = (1,2),  (2,1),  (2,2)  (276) 


In  general,  when  the  values  for  the  partials  in  Eqs  (274) 
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and  (275)  are  known,  the  following  equation  is  needed  to 

3Xmax(£v} 

obtain  one  element  of  »»  : 

3 F 1 

-y 


! 


3WSy>  5 " 2,  I 3W(V  3(”kl)ii  , . 

styab  ' £ it,  £ i?^a7  3fyab  ' 

3X »(F„) 


3y  performing  this  for  each  value  of  a and  b,  maXw  ^ ■ 

3 £y 

can  be  computed.  This  appears  to  be  a time  consuming  oper- 
ation, but  most  of  the  computation  is  iterative  and  can  be 
done  very  quickly  on  the  computer. 

3J(F  ) 

It  can  be  seen  from  Sq  (101)  that  — X-  is  required. 

3V 

Equation  (270)  gives  J(F ) for  the  multi-input  case  and,  as 
mentioned  earlier,  the  trace  operation  cannot  be  eliminated 
as  done  in  the  single  input  case.  Equation  (270)  becomes: 


N 


JCEy)  = Z ^ 


3a 


Ap~  ^ ( a ) B ( a ) 


-y— y 


_3 

3a  L 


Ap“*^(a)B(a) 


T T -1 
LC  R C 


(278) 


8 J(£ v) 

If  an  attempt  were  made  to  find  » directly,  eventually 

the  partial  of  a matrix  with  respect  to  another  matrix  would 
be  required.  To  avoid  this,  the  partial  is  evaluated  with 
respect  to  each  element  of  F . 

-y 
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3J(F  ) N m m m i 

3L-  = V jtr  2 XF„RF  iXiGiR“  C 

3 f . 3 f y y “ 

yab  yab 


+ XFRF  T — a XTCTR“1C  (279) 

(_**  y j y y gf 


where 


X = -2_[Ag-j(a)B(a) 
Da  L"F  “ _ 


(280) 


Equations  (272)  and  (273)  show  how  to  compute  — X.  Also 

Sfyab 


— F RF  T = YRF  T ♦ F RYT 

3 f -y — y — j -v — 

yab 


(281  ) 


where  Y is  an  mxr  matrix  consisting  of  all  zeros  except  for 

a one  in  the  a,b  position.  Performing  Eq  (279)  for  all  ele- 

3 J(F  ) 

ments  of  F gives  |£  . It  is  now  possible  to  compute 

3M  (FJ 


from  Eq  (101). 


m 0 i-v 

To  complete  the  algorithm  ij  is  required.  The  pro- 


cedure to  use  is  again  to  compute  first  • — *■■■  for  each  value 

DX  r. 


of  k.  From  Eq  (248) 


where  Eqs  (140)  and  (166)  can  be  used  to  compute 


Re- 


3 H 


3X  „ 

tv 


*• 


peating  Eq  (163) 


3 

“ -y 

3M(F) 
m — y 

T* 

T 

T 

— F = -Fc 

3 ipT 

^y  " 

(163) 


— Fc 


3 Mm  ( F„. ) 3F_ 

LX  and  - 

3X. 


Knowing  and  — ■■-X  , this  is  replaced  by 


’Ey 


3Mm(£y)  . £ £ 3 W 


3f. 


ab 


3 XT 


a=1  b=1  3 fy 


ax, 


(283) 


ab 


3M  (F  ) 

Performing  Eq  (283)  for  each  value  of  k then  gives  =~_jL.. 


This  is  the  required  gradient  to  maximize  M (F ) and  com- 

iu  jr 

pletes  the  modifications  needed  to  the  algorithm  equations 
of  Figure  5»  Table  XII  lists  the  differences  between  the 
single  and  multiple  control  algorithm. 

The  main  difference  between  the  single  input  case  and 
the  multi-input  case  is  the  additional  computational  burden 
resulting  from  multiple  inputs.  The  optimality  criterion 
and  general  approach  to  solving  the  problem  is  exactly  the 
same  for  both  cases.  It  can  be  seen  that  if  a problem  has 
multiple  unknown  parameters  as  well  as  multiple  inputs,  the 
computational  requirements  are  enormous  except  for  the  sim- 
plest cases. 
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v Conclusions  And  Suggestions  For  Further  Research 


In  this  Chapter  the  major  objectives  of  the  research 
as  well  as  the  contributions  are  discussed.  Also  dis- 
cussed are  areas  for  future  research,  all  of  which  are 
extensions  to  the  research  presented  in  this  disserta- 
tion. 

Objectives  And  Contributions 

The  objective  of  this  research  is  to  show  that  the 
estimation  of  unknown  parameters  in  the  plant  and  control 
matrices  of  a linear  system  can  be  enhanced  by  the  addition 
of  a constant  gain  output  feedback  control.  To  attain  this 
objective,  an  algorithm  is  developed  in  Chapter  II  that 
optimally  selected  the  feedback  vector  and  external  control 
sequence  for  a single  unknown  parameter,  single  input  con- 
trol problem.  An  energy  constraint  is  placed  on  the  external 
control  sequence  and  the  feedback  eigenvalues  are  confined 
to  a constraint  space.  This  space  is  selected  so  as  to  ful- 
fill stability  as  well  as  response  criteria  for  the  system. 

It  is  also  shown  that  for  the  case  in  which  the  number  of 
outputs  is  less  than  the  number  of  plant  states,  an  addi- 
tional constraint  is  placed  on  the  feedback  vector.  This 
additional  "consistency*'  constraint  is  required  in  order 
that  the  position  of  the  feedback  eigenvalues  could  be 
reached  by  feeding  back  the  outputs. 

The  optimization  criteria  that  are  chosen  are  raaximiz- 
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ing  the  trace  and  weighted  trace  of  the  Fisher  information 
matrix.  These  criteria  were  selected  because  they  increase 
the  information  about  the  unknown  parameters  in  the  output. 

By  properly  selecting  the  weighting  matrix,  maximizing  tr  SM 
is  equivalent  to  minimizing  the  determinant  of  the  inverse  of 
the  information  matrix.  Minimizing  the  determinant  reduces 
the  volume  within  contours  of  equal  probability. 

A closed  form  solution  to  the  optimization  problem  does 
not  exist,  so  an  iterative  technique  must  be  used.  The  tech- 
nique selected  is  a first  order  gradient  approach,  and  it  is 
used  to  find  the  values  of  £y0pt  and  ^ ^ because  of  the  guar- 
antee of  convergence  to  a local  solution  and  the  simplicity 
over  a higher  order  solution.  To  account  for  the  problem  of 
moving  the  eigenvalues  in  the  direction  that  maximized  the 
gradient  while  remaining  in  the  constraint  space,  the  gradi- 
ent projection  method  is  applied. 

In  Chapter  III  a single  unknown  parameter  example  is 
used  to  demonstrate  the  use  of  the  derived  algorithm.  The 
unknown  parameter  appears  nonlinearly  throughout  the  plant 
and  control  matrices.  A maximum  likelihood  estimator  is 


developed  and  used  to  estimate  the  unknown  parameter.  Sim- 
ulated system  outputs  are  corrupted  by  noise  from  a Gaussian 
noise  generator.  The  estimate  error  sample  means  and  stan- 


dard deviations  are  plotted  for  different  noise  levels  and 
control  steps.  The  plots  confirmed  that  the  addition  of 
feedback  controls  result  in  better  estimates,  especially  for 
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In  Chapter  IV  the  multiple  input  and  multiple  parameter 
cases  are  examined.  Both  result  in  substantially  increased 
complexity  of  the  problem,  and  the  computation  time  required 
to  solve  the  problem  could  become  a significant  limiting 
factor. 

It  is  difficult  to  select  a good  scalar  criterion  for 
the  multiple  parameter  cases.  The  more  appealing  c:  Lterion, 
because  of  the  simpler  computations  required,  could  result 
in  nonidentifiable  solutions.  The  criterion  that  would  pro- 
duce an  estimate  with  a minimum  error  covariance  usually  be- 
comes mathematically  intractable.  As  a compromise,  a com- 
bined criterion  was  selected.  This  ad  hoc  method  consisted 
of  using  the  maximization  of  the  trace  of  the  information 
matrix  as  the  criterion  as  long  as  the  determinant  of  the 
inverse  of  the  information  matrix  decreased.  If  this  did 
not  occur,  then  a weighted  trace  criterion  was  used,  where 
the  weighting  matrix  is  selected  such  that  the  determinant 
did  decrease.  This  required  that  the  feedback  matrix  remain 
constant. 

In  the  multiple  input  case  the  complexity  increased 
rapidly  as  the  dimension  of  the  input  controls  increased. 
Only  cases  in  which  the  dimension  of  the  input  controls  is 
small  would  be  practical  to  solve. 

Lopez-Toledo  (Ref  17)  addressed  the  feedback  control 
problem,  but  his  model  was  linear  in  the  parameters.  He 
also  did  not  address  the  problem  of  confining  the  feedback 
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system  eigenvalues.  A contribution  of  this  research  is  the 
solution  of  feedback  controls  in  which  the  parameters  are 
nonlinear  in  the  model.  Another  contribution  is  that  the 
solution  allows  constraints  to  be  placed  on  the  eigenvalues. 
The  computational  method  developed  in  this  research  reduced 
the  computer  time  required  to  solve  the  problem.  By  setting 
the  problem  up  in  a desired  format 

M = VT^V  * J(Fy)  (81+) 

maximizing  M can  be  accomplished  by  first  maximizing 


WVE 


J(F  ) 

-y 


i 

where  ^max^Ly)  is  maximum  eigenvalue  of  and  E is  the 
energy  constraint  on  V,  with  respect  to  F only.  Once  this 

y 

is  accomplished  then  the  controls  vector  V is  the  scaled 
eigenvector  corresponding  to  the  maximum  value  of  x _„(F  ). 

T+aX  "O' 

This  method  allows  the  feedback  matrix  to  be  solved  sepa- 
rately from  the  external  open-loop  controls. 

The  method  and  results  presented  in  this  research 
should  encourage  future  research  in  this  area.  Some  of  the 
areas  that  need  further  research  are  discussed  in  the  next 
section. 


Extension  For  Further  Research 

The  most  important  area  that  needs  further  research  is 
the  case  in  which  the  initial  conditions  of  the  states  are 


neither  zero  nor  free  to  be  chosen.  Since  the  research  re- 
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quired  zero  initial  states  or  freedom  to  select  them,  the 
system  would  have  to  be  initialized  to  these  conditions 
which  in  some  situations  may  be  neither  practical  nor  pos- 
sible. 

One  example  may  be  where  the  system  is  initialized 
properly  and  a certain  time  later  the  outputs  and  applied 
controls  are  used  to  estimate  the  parameters.  These  esti- 
mates could  then  be  used  as  a better  nominal  parameter  value 
for  future  runs,  but  the  system  is  probably  not  at  the  zero 
state.  Future  runs  could  not  occur  until  the  system  was 
reinitialized,  which  may  cause  a significant  delay. 

Another  area  for  research  would  be  to  investigate  mod- 
ern computational  techniques  to  decrease  the  computer  time 
required  to  obtain  a solution.  It  was  found  that  a signif- 
icant amount  of  time  was  required  even  for  the  simple  three 
state  problem.  In  order  to  apply  the  method  to  practical 
problems,  further  research  should  look  into  considering  the 
computational  requirements.  A rewarding  output  of  this  new 
research  may  be  the  application  of  this  technique  to  real- 
time problems. 

Another  area  would  be  to  investigate  the  possibility 
of  using  this  research  in  the  case  where  there  is  small  pro- 
cess noise.  Chapter  I mentioned  the  additional  complexity 
that  may  result  when  incorporating  process  noise  into  the 
system;  however,  further  research  may  show  that  for  small 
values  of  this  noise,  approximations  may  be  used  along  with 

the  procedures  developed  in  this  research. 
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One  area  that  was  briefly  considered  was  the  case  in 
which  the  value  of  the  measurement  noise  used  by  the  estima- 
tor was  different  from  the  real  value.  For  small  mlstuned 
R' s,  the  estimates  were  not  affected  substantially,  but  more 
time  was  required  by  the  estimator.  More  work  should  be 
spent  in  studlng  the  effects  of  mistuned  R' s since  this  ro- 
bustness issue  partially  answers  the  question  of  how  well 
this  method  will  work  on  real,  rather  than  simulated,  data. 

It  was  noticed  that  for  the  case  where  C-1  exists,  all 
the  eigenvalues  moved  to  the  boundary  of  the  constraint 
space  and  tried  to  cause  the  system  to  become  unstable. 

This  was  not  surprising  since  maximal  excitation  of  modes 
enhances  identification;  however,  attempts  at  proving  that 
thl3  in  fact  is  always  the  case  were  not  fruitful.  Further 
study  in  this  area  may  be  very  beneficial. 

The  final  suggested  area  for  further  research  is  to  de- 
termine if  there  is  an  optimal  set  of  columns  of  £ for  com- 
puting the  consistency  constraints.  There  may  be  a set  of 
columns  that  gives  constraints  that  produces  a better  input 
control.  Even  if' all  sets  yield  the  same  controls,  one  set 

may  minimize  the  computation  needed  to  compute  F . and 

— yopt 

2opf  Thi3  would  be  worth  investigating  further. 

The  purpose  of  the  dissertation  was  to  develop  an  op- 
timal feedback  control  that  will  enhance  the  Identification 
of  unknown  parameters.  This  was  accomplished  and  hope- 
fully it  will  Initiate  an  interest  in  applying  this  to  fu- 
ture problems. 
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Appendix  A 


Proof  That  £ tr 

j = l 


m rn  — I 

JaAP.C  j)hC*R  C 


Let 


G(  j ) = hiP/v(j)h 
E 


T -1 
CP  C 


( A-1  ) 
( A-2) 


Since  PA(j)  and  are  real,  symmetric,  positive  semidef- 
inite  matrices,  G and  E are  also  real,  symmetric,  positive 
semidefinite  matrices.  Let 


Y,  = tr  G(  j )E 

J 


(A-3) 


Since  G(j)  is  real  and  symmetric,  there  exists  an  orthogonal 
transformation  matrix  P such  that 


P G(j)P  = D(j) 

where  D(j)  is  a diagonal  matrix.  Now 

PTG(j)EP  = PTG(j)PPTEP 
= D(j)PTEP 
= D(j)J 


where  J = P^EP. 

For  all  x,  where  x is  a vector 


rn  <71  m 

x ~ Jx  = x* ( PxEP)x 


(A-4) 


(A-5) 
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= (Px)TE(Px) 


- 0 since  E * 0 


(A-6) 


which  implies  J ^ 0 . Since  the  trace  is  invariant  under 

a similar  transformation,  and  since  P is  orthogonal, 


tr  G( j )E  = tr  P“  G(j)EP 
= tr  PTG(j)EP 
= tr  D( j ) J 


* 2 dii3it 

i=l 


i i+1  n~ 

1 , 0 , • • • ,0  then 


(A-7) 


T r*1  i-1 

D(j),  J ^ 0 implies  d^,  jii  ^ 0 since  if  x = 0,---,0, 


0 ^ x Jx  = ji± 

0 * xTD(j)x  = di± 


(A-8) 

(A-9) 


tr  G(J)E  = 2 diiJ'ii  " 0 
i = 1 


(A-10) 


Therefore , 


yj  * 0 


(A-1 t ) 


n r rp  ti  _i  ~i  n 

2 tr  h PA(j)hC  R C = 2 Yi  - 0 


(A-1 2) 


This  holds  true  for  any  feedback  matrix  Fy. 


139 


Appendix  B 


Definition  of  a Vector  Derivative 

Let  the  scalar  x(^)  be  a function  of  the  n dimensional 
vector  By  convention,  the  definition  of  the  derivative 
is: 


This  is  sometimes  called  a gradient  and  physically  it  is 
finding  the  individual  gradients  of  x(^)  with  respect  to  a 
component  of  £ with  all  other  components  held  constant. 

In  the  n+1  space  Z,  where 


Eq  (3-1)  represents  a hyperplane  that  is  tangent  to  the  sur- 
face of  x(£)  at  the  point  where  the  derivative  was  evaluated. 


3 x(^)  f 3x(£)  a x(£) 


3x(£) 


( B— 1 ) 


Appendix  C 


Since 


3 a“F 


Ai(a)AF(a) 

3SL  F ”F 


" 3 ■«  _ ~]  _ j 3 A,-,(a ) 

= — -A^(a)  AF(a)  ♦ A^(a) — 

- 3 a J “ a a 


then 


Aj  + 1 = 

-A 


Ajj+1(a) 

3aJ+1(I) 


4+1(5) 


(C-4) 


(C-5) 


which  implies  that  it  is  true  for  the  case  k = j*1  , 

therefore,  it  is  true  for  all  k. 
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Appendix  D 


I 

1 / 


Proof  That  all  Available  Control  Energy  is  Utilized 

It  is  seen  from  Eqs  (22)  and  (83)  that  is  a real, 
symmetric,  positive  semidefinite  matrix.  Let  F be  constant 
and  let 

max  M = M1  = * J(F  ) (D-l  ) 

T 

Assume  0 < < E and  let  V-,  = bV^  where  b is  such 

that 

V^V2  = E (D-2) 

This  implies  that 

b^V,  . E (D-3) 

and 

b2  > 1 (D-4) 

Now 

m2  55  — 2— N— 2 + J%]  - M1  (D“5) 

From  the  fact  that  M1  = max  M 

M2  * YeHiSa  * 

■ »2v^v,  * j <Ey) 

> V^V,  * jfFy) 

> M,  (D-6) 

1 
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This  contradicts  Eq  (D-5)  so  the  assumption  that  0 < VjVj 
< E is  invalid. 

T 

Now  assume  V i V ^ =0  which  implies  that  = 0 . Eq 

(D-1)  now  becomes: 


max  M = M.  = J(F  ) 

J 


(D-7) 


Let  ^2-2  = ® and  since  Wjj  is  positive  semidefinite  there 
exists  a such  that 

M2  • vl**  * J(Fy) 

> J(V 

> M.  (D-8) 


This  also  contradicts  Eq  (D-5)  so  the  assumption  that 
T 

V1V1  =0  is  invalid.  Therefore, 
vfv,  = E 


(D-9) 


QED 
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(E-4) 


= BcklD  = BD^cklJ 


where  is  the  1-th  row  of  C. 
Therefore , 


— “Fa  «■  bftc~]d  = bdV 

3FL~  “ 

Since 


(E-5) 

(E-6) 


Appendix  F 


Derivation  of  Relationship  Between  F and  X 

-y  —■ r 

The  work  presented  in  this  appendix  was  developed  by 
B.  Porter  (Ref  25)  and  presented  by  Prof.  J.  D'Azzo  at  the 
Air  Force  Institute  of  Technology,  Wright-Patterson  AFB, 
Ohio  in  1978. 


Generalized  Control  Canonical  Form 

The  basic  structure  of  the  vector  state  equation  rep- 
resentation of  a system  can  be  incorporated  in  the  general- 
ized control  canonical  form.  This  canonical  form  exhibits 
inherent  properties  of  the  system  which  are  extremely  use- 
ful in  state  feedback  system  design.  This  form  is  readily 
amenable  to  closed-loop  system  design  which  achieves  a de- 
sired pole  placement.  While  it  is  not  the  only  form  which 
can  be  used  for  pole  placement,  it  has  advantages  for  high 
order  systems.  The  original  state  equation  for  the  open- 
loop  system  is 

x(kT*T)  = Fx(kT)  * Gu(kT)  (F-1) 


where  x(kT)  is  nxl , F is  nxn,  G is  nxm,  u(kT)  is  mxl , and 
the  rank  of  G is  m.  By  means  of  similarity  transformations 
of  x(kT)  and  u(kT)  and  by  the  addition  of  state  feedback, 
this  equation  can  be  converted  to  the  generalized  canonical 
form,  i.e.,  with  all  eigenvalues  located  at  the  origin,  of 
the  form 
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z(kT+T)  = F z(kT)  + G w (k'T) 

™ -C*-  “C  •'■c 


where 


Ic  = block  diagonal  (F^  , , • • •,  F^  ) 


m 


and 


G,  = block  diagonal  (g^  , g^.  ,• • , ^ 

12  ra 


a“-  ) 


(F-2) 

(F-3) 

(F-4) 


The  k^  in  equations  (F-3)  and  (F-4)  are  an  ordered  set  of 
integers  whose  sum  is  n and  are  known  as  the  control  invari- 
ants. These  control  invariants  are  uniquely  associated  with 
the  matrix  pair  (F,G)  and  are  identical  to  Kronecker's  min- 
imal column  indices  (kef  25)  for  the  singular  pencil  of  ma- 
trices (si  - F,G)  when  they  are  ordered  such  that 


and 


k1  - k2 


k^  * kg  + 


- K > 0 

m 


1 + k = n 
m 


(F-5) 

(F-6 ) 


The  matrices  F,  and  g,  are  of  size  k.xk.  and  k.xl  respec- 
i Ki  ill 

tively,  and  they  have  the  forms 


( 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

II 

•r| 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0_ 

1 

(F-7 ) 


As  shown,  the  submatrices  F,  are  comoanion  matrices  con- 

“Ki 
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I 


! i 
j,  | 


taining  ones  on  the  superdiagonal  and  the  remaining  elements 
are  zero,  and  the  vectors  are  null  except  for  a one  in 

its  k^  row  (which  is  its  last  row).  An  example  of  the  can- 
onical matrices  for  a system  which  has  k^  = 3 and  k^  = 2 
is 


0 

1 

0 1 

0 

0 

0 I 0 

0 

0 

1 1 

1 

0 

0 

o 1 0 

0 

0 

0 

0 

0 

G = 

1 ' 0 

— 

— \ 

— 

— c 

-4-  - 

0 

0 

0 1 

0 

1 

0 0 

1 

0 

0 

0 1 

0 

0 

1 

0 1 1 

(F-8) 


An  important  property  of  F is  that  it  has  an  index  of  nil- 
potency  (Ref  25)  cr  = 3,  which  corresponds  to  the  largest 
value  of  k^,  i.e.,  cr  = k^  =3* 

The  three  transformations  which  may  be  applied  to  the 
original  state  equation  (F-1 ) in  order  to  achieve  the  gener- 
alized canonical  forra(F-2)  are: 

(1)  A change  of  basis  in  the  state  variable  x(kT), 
that  is,  x(kT)  = Tz(kT),  where  det  T = 0.  Eq  (F-1) 
then  assumes  the  form 

z(kT+T)  = T*"^  FTz(kT)  ♦ T^GuOcT)  (F-9) 

(2)  A change  of  basis  in  the  control  variable  u(kT), 
that  is,  u(kT)  = Zv  (kT),  where  Z is  mxm  and  Det  Z = 0. 
Applying  this  transformation  changes  equation  (F-9)  to 
the  form 


1 
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f. 


z(kT«-T)  = T“'?Tz(kT)  ♦ T^GZv.  (kT) 
= T_1?Tz(kT)  - G v (kT) 


( F-1 0) 


(3)  The  introduction  of  feedback 


v (kT)  = w (kT)  - Lz(kT) 

“C  ^ C 


(F-1 1 ) 


where  L is  an  mxn  matrix  and  w,(kT)  is  an  mxl  input 
vector.  The  matrix  L is  required  to  convert  the  plant 
matrix  into  the  desired  form.  Substituting  v^(kT)  from 
Eq  (F-11)  into  Eq  (F-10)  yields  the  form 


z(kT  + T)  = T*"^  FTz_(kT)  + T“1GZwc(kT)  - T"1GZLz(kT) 
= (T~^FT  - T‘1GZL)z(kT)  + T^GZw  (kT) 

= F z(kT)  + G w (kT) 

“C  "* 


(F-1 2) 


It  may  not  be  necessary  to  apply  all  three  transformations 
in  every  case;  only  the  necessary  transformations  would  be 
used  to  achieve  the  canonical  form  of  the  state  equation. 
Also,  there  is  no  specified  order  in  applying  the  transfor- 
mations. It  should  be  noted  that  transformations  1 and  2 
produce  equivalent  matrices,  that  is,  the  control  invariants 
do  not  change  (Ref  25). 


Control  Invariants 

The  control  invariants  can  be  calculated  from  the  nx(nm) 
controllability  matrix 


g 3 [2.  £2.  f2g,  • * * i fH-’g] 


(F-1 3) 


It  has  been  3hown  by  Kalman  (Ref  14)  that  when  the  rank  of 
G is  m,  representing  a multi-input  system,  that  controlla- 
bility can  be  evaluated  from  the  nx(n-m+1)m  matrix 


2 = [G,  FG,  F2G, • • -,Fn"mG' 


The  columns  of  this  matrix  are 


(F-14) 


JT  — ["it  c*  . . . a-  F°"-  pJ  ...  Pit  F2*1.  . • • . F2-0" 

i L-1  ' a2’  »am>  £-»1  » E»2’  » £^m’  - aq  > » - am’ 

(F-1 5) 

For  a controllable  system  the  rank  ^ = rank  = n.  This 

means  that  the  matrices  £ and  £[  have  n linearly  independent 

column  vectors.  Define  the  column  vector  F1^,  to  be  regular 

J 

if  it  is  linearly  independent  of  the  columns  to  its  left  in 
A regular  basis  is  a basis  of  the  first  n regular  vec- 
tors contained  in  JJ.  These  are  rearranged  to  form 

2,  = {j»i  > > ' • * f £■]  > a2*  ' " — 2 E2 > ‘ • Ejh* 

• • F^m*"1^  (F-16) 

The  k^'s  are  the  control  invariants  of  the  system.  If  nec- 
essary, the  columns  of  G are  permuted  so  that  the  control 
invariants  are  in  decreasing  order,  as  given  by  Eq  (F-3). 

Transformation  to  Canonical  Form  - An  Example 


The  transformation  of  the  pair  (F,G)  to  the  canonical 

nair  (F  ,G  ) can  be  obtained  from  a digital  comnuter  using 
— c ’ — c 


f 


a coding  based  on  the  methods  such  as  those  of  Aplevich 
(Ref  2).  In  order  to  illustrate  the  transformation  using  a 
third  order  plant,  the  method  of  Luenberger  (Ref  18)  based 
on  his  second  canonical  form  is  illustrated  below  for  a sys- 
tem with  two  innuts: 


3 

5 

T 

~o 

T 

x(kT+T)  = 

0 

0 

i 

x(kT ) + 

0 

0 

0 

-2 

3 

1 

3 

u(kT) 


= Fx(kT)  + Gu(kT) 


x - 3 
0 
0 

Therefore 


-5  -1 

X -1 

2 X - 3 


= ( X - 1 )(  X - 2)(  X - 3) 


(F-17) 


(F-18) 


X-|=l»  X£  = 2,  x^  = 3 

In  reality  this  would  be  an  unstable  system,  but  as  an  ex- 
ample here  it  makes  no  difference. 


Thus  the  system  is  controllable  and  the  controllability 
index  (Ref  25)  u = 2.  The  first  three  columns  of  ^ are 

A 

linearly  independent,  so  ^ is  formed  as 


rank  ^ = rank  |jj, 

fg] 

E 

0 

1 

1 

r 

= rank 

0 

0 

, 

3 

= 3 

( F— 1 9 ) 

1 

II 

1 

3 

3 

9_ 

I — umm 
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a-  [a, 


*1  5 i -2J 


F -20 : 


1 3 I 3 


Forming  the  inverse  of  this  matrix  yields 


1 0 = e12 


1 -1  0 


( F— 21 


The  transformation  matrix  T-^  is  formed  by  using  the  bottom 

A-1 

row  of  each  partitioned  submatrix  of  ^ } as 


T”1  = e12F  = 0 0 1 

1 -1  0 


Now,  using  the  transformation 


(F-22 


x(kT)  = Tz(kT)  =100  z(kT) 


0 1 


(F-23 


The  plant  equation  becomes: 


z(kT+T)  = T“ ' FTz(kT)  ♦ T”'Gu(kT) 


= F z(k'T)  * G,u(kT) 

“C  C mm 


(F-24) 


In  order  to  complete  the  transformation  to  the  canonical 
form,  it  is  necessary  to  use  the  additional  transformations 

u(kT)  = Zv  (kT)  and  v (kT)  = w (kT)  - Lz(kT)  so  that  the 

— — —c  — c — c — 

state  equation  becomes: 


z(kT+T)  = (F  - GZL)z(kT)  ♦ (GZ)w  (kT) 

— — c — c — — —r.  — — r. 


— c-  -c 


= F z(kT)  + G w (kT) 

— c — — c — c 


(F-25) 


The  L matrix  is  of  order  2x3  and  the  Z matrix  is  2x2,  and 


hence  the  forms 


111  ■L'lZ  X13 


i21  122  123 


“ 11  “ 1 2 


“21  “22 


(F-26) 


Inserting  Z into  Eq  (F-24)  yields 


0 110 


z(kT+T)  = -2  3 1 0 z(kT ) 

8 0 | 3 


j 1 + ^“21  ^i  2 + 3 0)22  ^ ^ ^ F— 27 ) 


Since  G Z should  equal  the  generalized  canonical  form  G 
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r « 


— 

- 1 

— 

— 

0 

0 

0 

0 

1 

0 

• 

“11  * 3“21 

‘ u12  * *»22 

(F-28) 

0 

1 

_ “21 

“22 

Equating  coefficients  yields 


Z = 


1 -3 
0 1 


(F-29) 


The  matrix  Z will  always  be  in  the  upper  triangle  form.  In- 
serting L into 


z(kT+T)  = 


~0 

1 

~0 

~o 

0 

-2 

3 

0 

z(kT)  + 

1 

0 

8 

0 

3_ 

0 

1 

yields 


z(kT+T)  = 


w (kT)  - Lz(kT) 

— c — 

(F-30) 


0 

-2 

8 


1 

0 

0 

0 

0 

3 

0 

- 

1 

112 

b3 

0 

u 

121 

122 

X23_ 

= F z(kT)  ♦ G w (kT) 


z(kT) 


w (kT) 


(F-31) 


Since  F,  is  in  the  generalized  canonical  form 
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0 

0 

0 

0 

1 

0 

0 

1 

0 

xn 

X12 

X13 

= 

-2 

3 

0 

- 

0 

0 

0 

121 

122 

123_ 

8 

0 

3 

0 

0 

0 

Solving  for  L yields 


L = 


-2 

8 


3 0 

0 3 


The  three  transformations  used  are 


x(kT)  = 


1 0 1 
1 0 0 
0 1 0 


z(kT ) 


= Tz(kT) 


-3 

1 


u(kT)  = 

= Zv  (kT) 


v, (kT ) 

— G 


v (kT)  = w (kT) 

■"C  — c 


3 w, ( kT ) 


f-2  3 O' 

8 0 3, 

Lz(kT) 


z(kT) 


(F-3 2) 


(F-33) 


(F-34) 


(F-35) 


(F-36) 


which  transformed  the  original  system  representation  to  the 
generalized  control  canonical  form. 

Pole  Placement  by  State  Feedback 

Given  an  open-loop  plant  represented  by 
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x(kT+T)  = Fx(kT)  + Gu(kT) 


(F-1  ) 


the  intent  is  to  determine  a state  feedback  control  law 

u(kT ) = Kx(kT)  (F-37) 

which  will  produce  specified  eigenvalues  for  the  feedback 
system.  The  design  procedure  is  not  to  solve  for  K directly 
but  to  make  use  of  the  generalized  canonical  state  equation 
(F-11).  A feedback  represented  by 

w (kT)  = Hz(kT)  (F-38) 

produces 

z(kT+T)  = F z(kT)  + G Hz(kT) 

"""*  “"C  *""C  “ 

= (F  + G H)z(kT) 

= FjZCkT)  (F-39) 

where  F^  is  the  desired  feedback  system  matrix.  It  is  se- 
lected to  have  the  form 

F^  * block  diagonal  (c^  , ^ ^ ) (F-40) 

12  m 

where  each  c,,  is  a matrix  block  which  has  the  same  size  as 
“*i 

the  blocks  in  F and  is  in  companion  form.  In  Eq  (F-39)  the 

matrices  F , G , and  F,  are  now  known  and  H can  be  evaluated 

— c ’ — c —a  — 

by  equating  the  elements  of  (F  + G H)  and  F,.  A formal  ex- 

— c *— c—  —a 

pression  for  H can  be  obtained  since 


I 


This  leads  to 


H * S/d*  - FJ 


(F-42) 


In  order  to  return  to  the  original  state  space,  it  is  ob- 
served that 


u(kT)  = Zv^k'T) 

= Zw^kT)  - ZLz(kT) 
* Z(H  - L)z(kT) 

= Z(H  - L)T*"^  x(kT ) 

= Kx(kT) 


(F-43) 


Thus,  the  state  feedback  matrix  is 


K = Z(H  - L)T 


-1 


( F-44 ) 


The  advantages  of  using  the  generalized  canonical  form  are 
illustrated  by  continuing  the  example. 

Suppose  the  desired  feedback  system  eigenvalues  are 
A,  5 1 i *2  3 = *^±  order  to  ma^e  the  system  stable. 

The  factors  of  the  characteristic  polynomial  of  order  2 and 
1 respectively  are 


(x  - .5  + j.5)(  X-  .5  - j.5)  = X - X * .5 


(F-45) 


and 


X - 1 


(F-46) 


The  corresponding  matrix  blocks  of  the  desired  feedback  sys- 
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■ 

I 


. .a 


tem  matrix  F^  are  the  companion  form  matrices  associated 
with  these  characteristic  polynomials 


0 1 


1 -5  1 


s-k2  3 UJ 

From  Eq  (F-42) 


(F-47 ) 


( F-48 ) 


0 1 


0 0 1 


0 1 0 


-.5  1 0-0  0 0 
0 0 1 0 0 0 


1 0 


-.5  1 o 


-.5  1 o 
0 0 1 


(F-49) 


It  is  seen  that  the  rows  of  H are  made  up  of  the  coefficients 
of  the  characteristic  polynomials  and  zero  terms.  The  feed- 
back matrix  K is  then  determined  from  Eq  (F-44) 


K = Z(H  - L)T" 


1 -3|  (-.5  1 0 


-2  3 0 

8 0 3 


0 1 0 
0 0 1 


1 -1  0 
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Thus,  the  eigenvalues  are  1,  .5  ± j.5  as  desired.  This  con- 
cludes the  work  presented  on  the  generalized  canonical  con- 
trol form. 
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